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Section  I 


INTRODUCTION 


1.1  Program  Objectives 

The  Air  Force  Aircraft  Structural  Integrity  Program  (ASIP)  includes 
requirements  for  fracture  mechanics  analyses  of  aircraft  fleets  to  assess 
safety  limits  and  establish  inspection  intervals.  This  portion  of  the  ASIP 
includes  requirements  for  calculation  of  the  size  (as  a function  of  flight 
time)  of  cracks  which  might  be  present  in  fleet  airframes  [1] . The  initial 
conditions  for  the  calculations  are  cracks  of  prescribed  size  and  shape 
(specific  to  various  structural  details).  The  size  and  shape  values  are 
based  in  part  upon  observations  of  actual  airframe  cracks  from  past  service 
experience,  and  in  part  upon  the  degree  of  inspectability  of  the  airframe  [2]. 

The  crack-growth  calculations,  as  currently  required,  seek  to  assess  crack  size 
versus  time  for  the  deterministic  initial  conditions  outlined  above  and  for  load 
time-histories  also  assumed  to  be  deterministic.  The  load  spectra  are  prepared 
by  assembly  of  detailed  profiles  for  each  segment  (taxi,  climb,  cruise,  etc.) 
of  each  mission  type  (cargo  haul,  training,  air/air  combat,  etc.)  which  the 
fleet  is  expected  to  fly  [3]. 

In  reality,  variability  appears  in  the  actual  initial-crack  sizes,  in  the 
load  spectra  (through  individual  aircraft  usage  variation) , in  material  proper- 
ties, and  in  manufacturing  details.  Fleet  airframes  may  therefore  contain  a 
population  of  cracks  with  sizes  governed  by  a time- dependent  probability  dis- 
tribution. The  ASIP  criteria  seek  to  assess  the  effect  of  average  loads  on 
the  growth  of  above-average  initial  cracks.  Risk  analysis  is  an  associated 
structural  integrity  assessment  function  which  in  some  cases  seeks  to  quantify 
the  abnormal  possibilities  (extremely  large  initial  crack  size,  above-average 
loads)  in  terms  of  a probability  model.  Risk  analyses  of  the  latter  type  require 
a tremendous  volume  of  repetitive  calculations.  The  objectives  of  this  program 
were  to  assess  the  capabilities  of  hybrid  analog/digital  computers  to  perform 
such  analyses  in  an  efficient  manner,  and  to  employ  random-process  theory  to 
obtain  an  improved  description  of  the  random  components  of  loading  for  use  in 
risk  analysis. 

The  first  program  objective  was  organized  as  ten  specific  tasks,  nine 
involving  the  implementation,  verification  and  evaluation  of  an  analog  simulation 
of  one  or  more  aspects  of  currently  accepted  crack  growth-rate  models.  The 
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evaluation  included  an  assessment  of  the  ability  of  currently  availcible  Air 
Force  analog  computing  facilities  to  accommodate  the  simulations.  The  second 
program  objective  was  to  develop  a compatible  analog-based  method  for  simulating 
two-parameter  load  models  (i.e.  stress  peaks  and  ranges,  mean  and  alternating 
stresses,  or  stress  ranges  and  ratios)  based  on  power  spectra  and/or  peak- 
exceedance  curves  representative  of  an  aircraft  load  history. 

1.2  Modeling  the  Physical  Problem 

Formulation  of  a mathematical  model  of  the  physical  problem  of  crack  growth 
under  random  initial  conditions  and  random  loading  is  a difficult  and  extremely 
important  step.  Every  effort  must  be  made  to  insure  that  the  mathematical  model 
adequately  represents  the  physical  situation,  but  is  not  so  complicated  as  to 
preclude  numerical  calculations  because  of  excessive  costs.  Also,  the  model 
must  utilize  a data  base  which  is  very  large,  but  not  very  complete.  The  proce- 
dures used  to  develop  the  model  are  summarized  in  the  following  paragraphs. 

The  mathematical  model  can  be  separated  into  two  parts.  The  first  part  is 
the  environment  encountered  by  the  structure,  viz:  the  random  and  non-random 
components  of  loading.  The  second  part  is  the  crack-propagation  behavior  which 
results  from  the  loading. 

When  developing  the  loading  model,  every  effort  is  made  to  maintain  accuracy 
for  those  statistics  which  are  important  in  crack  propagation,  while  economies 
are  instituted  for  those  statistics  which  have  less  effect  on  crack  growth-rates. 
Numerous  investigators  in  the  past  have  treated  aircraft  loading  as  a Gaussian 
random  process.  In  a recent  investigation,  Shinozuka  [4]  considered  flight-by- 
flight loading  to  consist  of  ground  loads,  gust  loads,  GAG  cycles,  and  maneuver 
loads,  for  which: 


[The  experimental  evidence]  tends  to  indicate 
that  most  of  the  gust  loading  and  the  maneuver 
loading  for  cargo  or  transport  type  airplanes 
can  be  modeled  by  composite  Gaussian  processes 
while  the  maneuver  loading  for  fighter  type 
airplanes  can  be  approximated  more  closely  by 
single  Gaussian  processes." 


However,  other  evidence,  such  as  studies  of  air  turbulence  models  for  use  in 
flight  simulators  [5]  and  estimates  of  fighter  peak-n^-exceedance  curves  based 
in  part  on  flight  data  and  in  part  on  Vgh  calculations  [6]  provide  some  evidence 
of  non-Gaussian  characteristics  of  loading. 
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Crack-propagation  behavior  has  been  shown  to  depend  primarily  upon  the 
sequence  of  stress  minima  and  maxima  occurring  at  critical  locations  in  the 
structure,  coupled  with  the  instantaneous  size  of  the  fatigue  crack.  The  rate 
of  crack  growth  depends  primarily  on  the  difference,  As,  between  each  minimum 
and  its  succeeding  maximum.  The  Paris  equation,  the  simplest  form  of  crack- 
propagation  model,  shows  dependence  only  on  this  difference  [7].  This  is 
essentially  a linear  damage- summation  model  similar  to  the  Palmgren-Miner 
hypothesis  [8,9].  Linear  damage- summation  models  possess  an  appealing  simplic- 
ity especially  convenient  for  applications  to  random  loading,  in  that  the  sum 
of  crack-size  increments  for  many  load  cycles  is  independent  of  the  specific 
load  sequence.  Hence,  the  statistics  of  the  crack-increment  sum  can  be  related 
directly  to  the  time-average  statistics  of  a random  process  which  describes  the 
loading.  Forman  et  al.  [10]  have  included  a "layering"  effect  depending  on  the 
ratio  of  successive  values  of  stress  minima  and  maxima,  and  an  "acceleration" 
effect  depending  on  the  ratio  of  the  applied  crack-tip  stress  intensity  factor 
to  the  material  fracture  toughness,  providing  a greatly  improved  model.  With 
only  the  "layering"  effect,  this  is  still  essentially  a linear  damage- summation 
model.  With  the  incorporation  of  the  "acceleration"  effect,  the  time-average 
statistics  of  the  random  loading  are  no  longer  sufficient,  so  if  those  statistics 
are  nonstationary,  then  such  nonstationary  statistics  must  be  provided.  With 
this  required  statistical  description  of  the  loading,  most  of  the  appealing 
simplicity  of  the  linear  damage- summation  models  remains.  Simplified  random- 
process  models  for  damage  accumulation  have  been  developed  based  on  Miner’s 
rule  [11]  and  on  the  Forman  equation  for  crack-propagation  [12]. 

However,  the  work  of  Schjive  [13]  and  others  has  shown  that  fatigue  crack- 
ing in  structures  is  in  many  cases  sensitive  to  the  specific  load-application 
sequence.  This  effect  derives  from  the  material  "memory"  associated  with  plastic 
yielding  and  subsequent  residual  stresses  near  the  crack  tip.  Modifications  of 
the  linear  damage- summation  models  for  crack-propagation  have  been  made  by 
Wheeler  [14] , Willenborg  et  al  [15] , Dill  and  Saff  [16]  , and  others  to  account 
for  material  memory  effects.  Past  implementations  of  these  models  have  generally 
required  cycle-by-cycle  summation  of  crack-size  increments.  Several  efficient 
digital  computer  codes  (e.g..  Ref.  17)  have  been  developed  for  this  purpose. 

More  recently,  Galagher  [18]  has  shown  for  transport  aircraft  that  flight-by- 
flight  summation  can  be  used  to  predict  crack  growth  accurately  if  the  rate 


3 


equation  is  formulated  on  a per-f light  basis,  with  parameters  fit  to  flight-by- 
flight  crack-propagation  test  data. 

With  regard  to  load  models,  it  was  the  judgment  of  the  present  investigators 
that  load  models  should  not  be  restricted  to  Gaussian  load  spectra  unless  both 
adequate  modeling  accuracy  and  simplification  of  the  resulting  analysis  could  be 
assured.  As  will  be  described  in  Subsection  2.2,  such  assurance  is  lacking.  On 
the  other  hand,  the  crack-propagation  models  described  above  have  been  repeatedly 
verified  by  comparison  of  predictions  with  laboratory  experiments  using  determin- 
istic loading.  Hence,  the  present  investigators  accepted  these  crack-propagation 
models  in  this  program. 

1.3  Results  of  Previous  and  Present  Programs 

The  statistics  of  crack  size  were  investigated  by  Monte  Carlo  simulation, 
using  a small  analog  computer,  in  a previous  program  under  Air  Force  Contract 
F33615-74-C-3046.  Repeated  runs  with  the  loading  simulated  by  a white-noise 
generator  filtered  through  analog  sample-and-hold  circuitry  yielded  crack-size 
histograms  to  which  3-parameter  Weibull  distributions  were  fit  by  maximum- 
likelihood  estimation  [19].  Crude  models  were  used  for  both  the  crack  growth- 
rate  behavior  and  the  random  loading  because  the  primary  objective  of  the 
investigation  was  to  assess  the  feasibility  of  performing  Monte  Carlo  simula- 
tion on  an  analog  computer  with  the  rising  exponential  type  of  dynamics  generally 
exhibited  by  crack  size  versus  time.  A simplified  form  of  the  Paris  equation 

without  local  geometry  effects  was  taken  as  the  growth-rate  model: 

N 

da/dN  C (ASi/a)  ^ (D 

P 


where 

a = Instantaneous  crack  size 

N = Number  of  applied  stress  cycles 

As  = s - s . = Stress  range 

max  min 


and  where  C^,  are  empirical  constants.  For  the  purpose  of  the  simulation, 

was  assumed  equal  to  4 (representative  of  a wide  range  of  aluminum  alloys) 

and  C was  fit  to  experimental  data  for  a 7075-T6  alloy  [20].  The  stress  range, 
P 

As,  was  assumed  as  a Gaussian  random  process  with  parameters  fit  approximately 
to  available  fighter  exceedance  data  [6] . The  investigation  did  demonstrate 
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the  feasibility  of  analog  Monte  Carlo  simulation  for  cracks  grown  by  random 
loading.  However,  the  specific  results  obtained  were  later  found  to  be  in 
error  due  to  an  unexpected  low-frequency  resonance  in  the  sample-and-hold 
circuitry  which  was  discovered  in  the  course  of  the  present  program. 

In  the  present  investigation,  eight  of  nine  specific  simulation  tasks  were 
developed  as  analog  programs  which  were  implemented  on  a CI-5000  hybrid  analog/ 
digital  computer  at  the  USAF  Aeronautical  Systems  Division  (ASD/ADSD) . Each 
task  was  demonstrated  by  means  of  a few  verification  runs.  However,  full 
Monte  Carlo  simulations  were  not  performed  for  reasons  which  arose  from 
developments  under  the  remaining  task  in  the  investigation. 

Development  of  this  latter  task  (two-parameter  simulation  of  random 
loads)  led  to  the  discovery  of  the  low-frequency  sample-and-hold  resonance 
mentioned  above.  The  use  of  sample-and-hold  circuits  was  therefore  rejected. 

Also,  considerations  of  the  limits  of  analog  equipment  for  processing  high- 
frequency  signals  led  to  the  conclusion  that  much  more  computation  time  would 
be  required  for  one  airplane  life  than  had  previously  been  estimated.  Hence, 
the  analog  Monte  Carlo  simulation  approach  to  risk  analysis  appeared  to  be  much 
less  efficient  than  previously  expected. 

In  the  meantime,  other  developments  under  the  load  modeling  task  led  to 
two  ideas  for  improvement  of  the  simulations.  First,  the  crack  growth-rate 
equations  were  reformulated  in  terms  of  a damage  parameter  which  varies  inversely 
with  crack  size.  This  transforms  the  dynamics  of  crack  growth  from  exponential 
to  nearly  linear  behavior.  The  new  formulation  was  used  to  implement  the  eight 
specific  analog  simulation  tasks. 

Second,  a general  study  of  possible  approaches  to  the  simulation  of  random 
loads  led  to  the  concept  of  applying  estimation  theory  to  the  crack-growth  predic- 
tion problem.  Estimation  theory  was  developed  by  modern  theorists  in  the  field 
of  guidance  and  control,  primarily  as  a tool  for  optimizing  designs  for  electronic 
control  systems  [21,22].  The  present  study  has  indicated  that  estimation  theory 
can  be  applied  to  the  crack-growth  prediction  problem  to  obtain  the  sought-after 
computational  efficiency  for  risk  analysis.  Hence,  the  application  of  estimation 
theory  to  crack-growth  dynamics  was  considered  as  a part  of  the  random- load  simula- 
tion task. 

Details  of  the  developments  made  under  this  program  are  presented  in  the 
following  sections.  Section  II  discusses  the  transformation  from  crack  size  to 
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damage  parameter,  and  the  benefits  obtained  by  using  a damage-parameter  formula- 
tion for  risk  analysis  of  random  loads  and  crack  sizes.  In  Section  III,  the 
implementations  of  the  specific  analog  simulation  tasks  are  presented,  together 
with  example  runs  and  evaluations.  Section  IV  presents  conclusions  and  recom- 
mendations, including  discussion  of  the  application  of  estimation  theory  to 
the  crack-growth  prediction  problem. 
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Section  II 


FORMULATION  FOR  COMPUTATION  OF 
CRACK-GROWTH  STATISTICS 


2.1  The  Damage  Parameter  ”h" 

For  analysis  of  statistics  of  the  crack  size  resulting  from  the  application 
of  random  loading  to  a fatigue  crack,  it  can  be  very  advantageous  to  use  a damage 
parameter  "h”  in  place  of  the  crack  size  "a".  For  some  mathematical  models  of 
crack  propagation  behavior,  such  as  the  Paris  equation  with  no  geometry  factor, 
an  appropriately  chosen  definition  of  "h"  can  simplify  the  statistical  deriva- 
tion to  triviality.  For  other  more  complicated  models,  the  statistical  deriva- 
tion is  at  least  somewhat  simplified. 

The  simpliest  commonly  used  model  of  crack  propagation  behavior  is  the 

Paris  equation  (see  Eq.  1) . If  the  growth- rate  exponent  N is  chosen  equal  to 

P 

4 (representative  of  aluminum  alloys) , then: 

da/dN  ^ C^(AS)^  a^  (2) 

If  the  damage  parameter  is  simply  chosen  as 

h = 1/a  (3) 


then  substitution  of  Eq.  3 into  Eq.  2 gives: 

dh/dN  ::  -c  (AS)^  (4) 

P 


The  change  in  h then  depends  only  on  the  applied  loading,  and  not  on  the  value 
of  h itself.  For  the  more  general  case  of  N^^2 , the  damage  parameter  can  be 
chosen  as 


h 


1-N  /2 
P 


(5) 


which  provides 


dN 


(1 


N 

C (AS) 

2 P 


N 

P 


Cases  with  N <2  are  considered  later  in  this  subsection. 
P" 


(6) 
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Figure  1 shows  constant-stress-amplitude  time-histories  of  crack  size  "a" 

and  damage  parameter  "h",  for  crack  growth  according  to  the  Paris  model.  These 

histories  show  some  problems  which  will  affect  analog  simulation  using  crack 

size.  Simulation  of  crack  size  will  always  saturate  the  analog  equipment  short 

of  complete  structural  failure,  as  shown  at  (A) . Scaling  crack  size  to  have  its 

critical  value  a below  saturation  level  (B)  makes  the  initial  size  quite  small 
cr 

(C) , and  hence  subject  to  analog  error.  This  error  is  important  because  da/dN 
is  small  initially,  so  that  a very  small  error  in  the  voltage  representing  initial- 
crack  size  will  yield  a very  large  error  in  estimated  time  to  the  critical  condi- 
tion. This  time-history  is  also  subject  to  analog  error  at  later  times,  because 
of  the  extreme  rate  of  change  da/dN,  but  this  form  of  error  will  have  much  less 
effect  on  the  estimate  of  time  to  the  critical  condition. 

For  the  simple  Paris  model,  an  appropriate  definition  of  the  damage  parameter 
**h'*  has  eliminated  the  dependence  of  the  rate  of  damage  accumulation  on  the 
instantaneous  value  of  accumulated  damage.  However,  the  advantages  of  the  use 
of  a damage  parameter  are  reduced  for  more  complicated  and  more  realistic  models 
of  crack  growth  behavior,  such  as  the  Forman  equation. 

The  Forman  equation  [10]  is  usually  written  as: 


da 

dN 


C (AK) 
F 

(l-R)K  - 
c 


N 


F 


AK 


(7) 


where 


AK  - As/a 


(8) 


if  local  geometry  effects  are  neglected.  The  parameters  C^,  N^,  and  are 
material  constants.  The  crack  size,  a,  and  the  stress  difference  As  - "’min^ 

are  the  same  as  for  the  Paris  equation.  The  stress  ratio  R = i^^^^duces 

the  effect  of  the  individual  values  s^.  and  s . 

mxn  niaA 

It  is  convenient  to  substitute  for  AK,  cind  to  use  As  and  as  the  describers 

of  the  stress.  This  provides  the  following  equivalent  form  of  the  Forman  equation: 


(AS) 


N^/2 

a 


AS 


-]  [ 


s -s 

c max 


(9) 
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where 


S = K //a 
c c 


(10) 


Equation  9 may  be  considered  in  three  parts.  The  first  is  the  basic  relation 


da 

dN 


K 

c 


% 

(ASi/a) 


(11) 


which  is  like  the  Paris  equation.  The  second  and  third  parts  are  multiplicative 

corrections  for  experimentally  observed  effects:  the  "layering"  effect  (s  /As) , 

max 

cind  the  "acceleration"  effect  [S  /(S  -s  )]. 

c c max 

A damage  parameter  may  be  defined  for  the  Forman  equation  in  the  saime  manner 

as  was  done  for  the  Paris  equation.  Assuming  N >2, 

F 


l-Np/2 

h = a 


(12) 


dN 


S 

max 
^ AS 


"S  -S 
c max 


(13) 


where 

l/(Nj,-2) 

" V (14) 

This  particular  definition  of  the  damage  parameter  results  in  a great  simplifica- 
tion of  the  crack  growth  relationship,  but  fails  to  achieve  the  triviality  found 

for  the  Paris  model.  The  acceleration  effect  [S  / (S  -s  ) ] still  depends  on  the 

c c max 

instantaneous  value  of  the  damage  parameter.  However,  this  dependence  is  fairly 
simple  and  well-behaved.  Analysis  of  the  statistics  of  final  crack  size  should 
still  be  manageable. 

Figure  2 shows  constant-stress-amplitude  time-histories  of  crack  size  "a" 
and  damage  parameter  "h",  for  crack  growth  according  to  the  Forman  model.  The 
time-history  of  crack  size  possesses  the  same  problems  as  for  the  Paris  model. 

The  time  history  of  the  damage  parameter  exhibits  one  problem  not  seen  with  the 
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Paris  model:  large  rate  dh/dN  near  complete  failure.  However,  this  problem 

is  not  very  important;  it  compares  with  the  least  of  the  problems  cited  for 

crack-size  history,  and  is  less  severe.  It  should  also  be  noted  that  if  h^^ 

is  based  on  an  isolated  high  overload,  the  rate  dh/dN  will  not  increase 

excessively  until  well  after  h is  exceeded.  This  is  because  h so 

chosen  provides  almost  zero  probability  of  exceedence  of  (catastrophic 

structural  failure)  for  h > h (aircraft  still  in  service) . As  a result, 

cr 

there  will  be  small  probability  of  s near  S , and  consequently  the  crack- 

max  c 

growth  acceleration  will  not  be  severe. 

Although  the  above  choice  of  the  damage  parameter  is  extremely  useful,  a 
definition  of  a damage  parameter  which  would  eliminate  all  dependence  of  crack 
growth-rate  on  current  damage  would  be  especially  desirable.  Unfortunately, 
incorporation  of  the  acceleration  effect,  which  includes  the  value  of  crack 
size  in  a very  different  fashion  from  the  first  two  effects,  makes  this 
impossible.  However,  there  may  be  some  methods  of  further  simplifying  the 
statistical  analysis  by  changing  the  format  of  the  crack-growth  model.  It  is 
not  necessary  to  use  the  Forman  model,  if  another  relationship  can  be  found 
which  provides  similar  agreement  with  the  experimental  data.  For  example, 
acceleration  need  not  be  modeled  as  a multiplicative  effect.  It  may  be 
additive,  in  which  case  the  crack  growth  rate  may  be  modeled  as: 


da  _ ^F 
dN  K 

c 


(AS,^)  ^ 


'■  As 


s 

c max 


(15) 


If  the  same  damage  parameter  "h"  as  described  for  the  Forman  equation  is  used, 
the  results  are  the  same  except  that  the  term  depending  on  the  current  damage 
is  additive  rather  than  multiplicative.  Such  a chainge  would  further  simplify 
the  statistical  analysis  of  final  crack  size.  The  use  of  Eq.  15  has  not  been 
pursued  in  this  study  because  additive  acceleration  models  have  not  been 
verified  by  comparison  with  experimental  results. 

The  definitions  of  the  damage  parameter  "h"  presented  earlier  in  this 
subsection  depended  on  values  of  or  greater  than  two.  For  some  materials, 
such  as  steel,  these  values  may  be  equal  to  or  less  than  two,  requiring  different 
definitions  of  h.  These  definitions  will  be  considered  here  using  the  Paris 
equation,  and  may  be  easily  extended  for  application  to  the  Forman  equation. 
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When  Np  is  less  than  two,  the  mathematical  definition  of  the  damage  param 
eter  is  identical  to  that  for  Np  greater  than  two  (Eg.  5) . However,  the  behav 
ior  of  the  damage  parameter  will  be  different,  in  that  it  will  increase  rather 
than  decrease  from  its  initial  value,  as  shown  in  Fig.  3.  When  Np<2,  "h"  will 
reach  infinity  when  crack  size  "a”  reaches  infinity,  so  it  is  clear  that  for 
the  ideal  case  of  a semi-infinite  plate  complete  structural  failure  will  not 
occur  in  finite  time.  This  compares  with  the  earlier  case  of  Np>2,  where  "h" 
reached  zero  when  "a"  reached  infinity,  and  complete  structural  failure  did 
occur  in  finite  time  for  the  ideal  case. 

When  Np  is  equal  to  2 , Eg.  5 is  singular,  and  an  entirely  new  definition 
must  be  found.  If  the  damage  parameter  is  defined  as 

h - Inia)  (15) 


then  the  Paris  equation  becomes 


dh 

dN 


C (As) ^ 
P 


(17) 


As  in  the  case  of  Np<2,  h increases  with  time,  and  an  infinite  value  of  h 
represents  complete  structural  failure. 

No  matter  what  the  value  of  Np,  it  is  not  necessary  to  utilize  the  defini 
tion  of  "h"  which  reduces  the  Paris  equation  to  an  exactly  constant  rate  of 
change  of  "h”.  Such  a definition,  as  presented  for  various  values  of  Np  so 
far  in  this  subsection  will  be  referred  to  as  the  "natural"  definition  and  is 
desirable  in  that  it  simplifies  the  ensuing  analysis.  However,  if  the  defini- 
tion of  the  damage  parameter  is  chosen  for  a value  of  Np  different  from  but 
close  to  the  actual  value,  the  analysis  will  remain  sufficiently  simple.  For 
the  case  of  the  Forman  equation,  where  the  rate  of  change  of  h is  not  exactly 
constant  even  for  the  natural  definition  of  h,  there  may  be  no  noticeable 
increase  in  analytical  difficulty.  In  fact,  the  analysis  with  the  Forman 
equation  may  even  be  simplified  by  defining  the  damage  parameter  so  that  the 
time  rate  of  change  of  h is  as  constant  as  possible  including  consideration 
of  the  layering  and  acceleration  effects.  This  is  accomplished  by  defining 

h as  if  N_  were  somewhat  larger  than  its  actual  value.  It  cannot  be  over- 
r 

emphasized  that  the  definition  of  h utilized  for  a particular  analysis  is  a 
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free  choice  by  the  analyst,  and  its  relation  to  the  natural  definition  of  h 
for  the  particular  crack  growth  rate  equation  is  guided  only  by  the  analyst’s 
desire  to  obtain  a simple  and  accurate  analysis. 

One  reason  for  utilizing  definitions  of  h different  from  the  natural  defi- 
nitions, even  for  the  Paris  equation,  is  for  comparison  of  behavior  for  various 
values  of  N^.  Figure  4 shows  time  histories  for  the  damage  parameter  defined 
as  h = In (a),  the  natural  definition  for  = 2.  Time  histories  are  shown  for 
behavior  with  values  of  Np  equal  to  1.8,  2.0,  and  2.2.  The  three  cases  shown 
have  identical  initial  crack  size,  and  their  Paris  equation  coefficients  Cp 
were  chosen  to  provide  identical  crack  growth  rate  at  the  initial  time. 

2.2  The  Statistics  of  Random  Loading 

As  described  in  Subsection  1.2,  the  loading  sequence  characteristics  which 

have  the  greatest  effect  on  crack  growth  are  the  values  of  minima  and  immediately 

succeeding  maxima  on  the  applied  stress,  s . and  s . The  most  important  effect 

min  max 

is  that  of  As  = s - s . , with  the  individual  values  having  a secondary  effect, 

max  min 

The  statistical  model  of  the  environment  encountered  by  the  structure  should  there- 
fore emphasize  these  characteristics,  even  though  the  corresponding  statistics 
might  be  much  more  difficult  to  establish  than  other  characteristics  of  the  loading. 

The  Paris  model  of  crack  growth  utilizes  only  the  values  of  As.  For  this 
model,  a probability  density  function  of  this  one  variable  suffices  as  a complete 
loading  model.  The  Forman  model  and  other  reasonably  realistic  models  include 

the  effect  of  the  individual  values  of  s . and  s on  the  crack  growth  rate. 

min  max 

Because  of  the  form  of  the  crack-growth  models,  individual  probability  density 

functions  of  these  two  variables  are  not  sufficient.  At  least  some  information 

about  their  correlation  is  necessary,  and  a joint  density  function  of  the  two 

variables  would  be  especially  desirable. 

There  are  three  approaches  which  may  be  used  to  determine  the  desired 

statistical  information  about  s . and  s . The  first  approach  assumes  that 

min  max 

the  loading  is  Gaussian  noise,  and  attempts  to  analytically  derive  the  desired 

statistical  information.  Given  Gaussicin  loading,  s,  the  individual  probability 

density  functions  of  s . and  s are  well  known  [23] . However,  the  probability 

min  max 

density  of  As  = s - s . depends  upon  the  correlation  between  each  minimum 
max  min 

and  its  associated  (immediately  succeeding)  maximum,  and  can  be  obtained  only 
with  difficulty  by  a numerical  integration  procedure  developed  by  J.R.  Rice 
[24,25,26].  It  should  be  noted  that  As  by  definition  cannot  be  negative,  and 
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a normal  probability  density  always  provides  nonzero  probability  of  negative 

values.  Hence,  even  when  the  loading  is  Gaussian  noise,  the  probability  density 

function  of  As  cannot  be  normal.  An  example  of  the  probability  density  function 

of  As  is  shown  in  Fig.  5,  for  a particular  Gaussian  noise  loading. 

The  joint  density  function  of  s . and  s , desired  for  more  realistic 

min  max 

crack  growth  models,  is  much  more  difficult  to  find.  The  necessary  procedures 

for  developing  such  a joint  density  function  are  available  [24].  Unfortunately, 

the  procedures  include  numerical  methods  because  closed-form  expressions  could 

not  be  formulated.  Utilization  of  the  procedures  requires  a substantial  effort 

which  must  be  repeated  for  each  new  form  of  power  spectral  density  considered. 

For  the  present  analog  simulation  effort,  maximum  use  is  made  of  Rice's 

results  without  utilization  of  these  procedures.  An  interim  solution  for 

modeling  time-histories  of  s . and  s has  been  developed  for  use  in  an  analog 

min  max 

Monte  Carlo  approach.  Although  many  approximations  and  assumptions  are  made, 
for  simplicity  and  because  of  the  scarcity  of  information,  the  interim  model 
does  produce  time-histories  whose  statistics  approximate  the  density  functions 
shown  in  Fig.  5 reasonably  well. 

The  second  approach  to  the  determination  of  the  statistics  of  s and 

min 

Smax  utilizes  a Monte  Carlo  analysis.  Gaussian  noise  with  a particular  power 

spectral  density  is  produced  on  analog  computer  equipment  by  driving  a linear 

system  with  a white-noise  generator.  The  digital  portion  of  a hybrid  analog/ 

digital  computer  then  samples  the  Gaussian  noise  sufficiently  rapidly  to 

determine  the  occurrences  and  values  of  minima  and  maxima.  These  values  are 

processed  by  the  digital  computer,  over  a substantial  time  period,  to  determine 

the  s . and  s statistics.  This  may  not  only  be  more  cost-effective  than 
min  max 

Rice's  derivation  for  the  straightforward  Gaussian  case,  but  it  is  much  more 
suitable  for  extension  to  more  complicated  cases.  For  instance,  certain  forms 
of  non-Gaussian  random  loading  may  be  treated,  if  the  loading  characteristics 
are  well  enough  known  for  appropriate  signals  to  be  generated  for  sampling  and 
statistical  processing. 

The  third  approach  to  the  determination  of  the  statistics  of  s . and  s 

min  max 

is  the  most  accurate,  but  depends  on  the  availability  of  extensive  actual  load 
histories.  This  approach  involves  examination  of  the  load  histories  to  deter- 
mine the  successive  values  of  s . and  s , and  inspection  of  these  values  to 

min  max 

determine  their  statistics.  The  loading  is  not  necessarily  Gaussian  noise;  in 
fact,  the  investigator  need  not  have  any  a priori  knowledge  of  the  loading 
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statistics.  This  makes  this  third  approach  the^most  flexible  of  all.  The  steps 
involved  are  represented  in  Fig.  6.  As  shown,  available  flight  data  will  give 
sampled  values  of  flight  parameters  at  discrete  times.  These  data  can,  for 
example,  be  interpolated  to  provide  assumed  time  histories.  The  maxima  and 
minima  can  then  be  identified,  and  the  important  characteristics  such  as  As 
can  be  established  for  each  loading  cycle.  (A  simpler  but  still  useful  approach 
is  to  program  a digital  simulation  of  a discrimination  counter.)  These  results 
can  then  be  inspected  to  determine  the  statistics  of  the  input  to  the  crack- 
growth  dynamic  system. 

In  previous  investigations  of  the  statistics  of  aircraft  loading,  there 

Tiave  been  attempts  to  characterize  the  loading  sequence  in  some  analytically 

describable  fashion  [4,5].  The  stress  "s"  has  been  assumed  to  be  a Gaussian 

processes.  This  might  appear  to  be  a simplification;  i.e.,  an  approximation 

which  makes  the  subsequent  utilization  of  the  load  model  easier.  However, 

there  is  no  resulting  simplification  because,  as  described  above,  it  is  extremely 

difficult  to  derive  the  statistics  of  s . , s , etc.  from  the  power  spectral 

min  max 

density  of  s,  and  the  resulting  statistics  are  not  Gaussian  anyway.  Therefore, 

the  third  approach  described  above  for  determination  of  the  statistics  of  s . 

min 

and  s (i.e.,  direct  inspection  of  large  quantities  of  flight  data)  is  the 
max 

method  of  choice  for  crack  growth  risk  analyses.  This  approach  avoids  both  the 
necessity  to  analytically  describe  the  loading,  and  the  difficult  analytic 
determination  of  its  important  statistics. 

The  crack-growth  dynamic  system  can  be  modeled  with  either  load  cycles  or 
time  (flight  hours)  as  the  independent  variable.  If  the  system  is  modeled  with 
load  cycles  as  the  independent  variable,  some  correspondence  must  be  established 
between  time  and  cycles.  Fortunately,  the  number  of  load  cycles  per  unit  time 
varies  little  among  various  samples  of  the  random  loading,  even  though  the  impor- 
tant characteristics  of  the  loading  may  vary  substantially.  Sufficient  accuracy 
might  therefore  be  obtained  by  utilizing  the  mean  value  of  load  cycles  per  unit 
time  to  transform  from  the  statistics  in  terms  of  load  cycles  produced  by  the 
dynamic  model,  to  statistics  in  terms  of  flight  hours.  However,  when  direct 
inspection  of  flight  data  is  used  to  obtain  loading  statistics,  there  is  no 
need  to  depend  upon  the  expected  number  of  load  cycles  per  unit  time.  Since 
the  sampled  flight  data  includes  time  information,  the  statistics  can  simply 
be  gathered  in  terms  of  occurrences  per  unit  time.  Expected  total  number  of 
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load  cycles  per  unit  time  may  be  one  of  the  statistical  parameters  measured, 
but  it  need  not  be  a basis  on  which  all  other  statistics  are  dependent  for 
their  expression  in  terms  of  time. 

Consideration  of  the  statistics  necessary  for  accurate  determination  of 
the  "acceleration  effect"  provides  further  reason  for  obtaining  statistics 
from  direct  inspection  of  flight  data.  The  Forman  equation  (Eq.  13)  indicates 
what  the  important  characteristics  of  the  loading  will  be  for  the  basic,  layer- 
ing, and  acceleration  effects.  The  basic  and  layering  effects  merely  require 

some  moments  of  the  joint  probability  density  function  of  s . and  s , or 

min  max 

perhaps  that  of  As  and  s for  better  accuracy  of  results.  The  acceleration 

max 

effect  will  require  adequate  information  concerning  the  occurrences  of  excep- 
tionally large  values  of  s .It  is  known  that  a simple  Gaussian  model  for 

max 

stress  "s"  will  provide  inaccurate  estimates  of  occurrences  of  such  large  values 

of  s [5] . Improved  analytically  describable  models  would  provide  better 
max 

estimates,  but  their  accuracy  would  always  be  questionable.  Direct  inspection 
of  flight  data  can  avoid  this  problem,  as  long  as  the  need  for  the  statistics 
of  large  values  of  is  predetermined.  The  first  few  moments  of  the  prob- 
ability density  of  cannot  be  expected  to  provide  a sufficiently  accurate 

prediction  of  the  occurrences  of  exceptionally  large  values,  but  separate 
information  concerning  level  exceedances  can  be  gathered  during  inspection  of 
the  flight  data.  If  the  acceleration  effect  is  multiplicative  as  in  the  Forman 

equation,  then  joint  statistics  with  As  and  s must  be  considered  for  the 

max 

level  exceedances.  However,  if  the  acceleration  effect  can  be  modeled  as 
additive,  then  the  statistics  of  level  exceedances  may  be  considered  separately 
from  the  statistics  of  As  and  greatly  reducing  the  complexity  of  the 

statistics  which  must  be  handled.  In  order  to  determine  whether  an  additive 
acceleration  effect  models  the  material  behavior  sufficiently  well,  empirical 
data  such  as  presented  in  Ref.  20  should  be  consulted. 

Interaction  effects  between  load  cycles  must  also  be  considered.  Load- 
interaction  consists  primarily  of  a reduction  in  crack  growth-rate  during  the 
few  cycles  immediately  following  the  occurrence  of  a very  large  maximum  stress. 
When  performing  cycle-by-cycle  analyses,  this  effect  can  be  modeled  by  various 
retardation  rules,  of  which  the  most  accurate  [15,16]  require  current  crack-size 
values  for  the  analysis.  However,  the  effects  of  retardation  can  still  be 
incorporated  in  the  load  model  if  suitable  assumptions  about  crack  size  are 
introduced. 
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An  illustrative  example  based  on  the  Willenborg  model  will  be  given  here. 

The  Willenborg  model  [15]  replaces  the  applied  crack-tip  stress  intensity  factor, 
K,  by  an  effective  factor  defined  as: 

K*  = K-K  ' (18) 

R 


where 


.OL 


K = K'"  A - Aa/z 
R max  '-’ll 


(19) 


and  where  is  the  maximum  stress  intensity  factor  associated  with  a prior 

max 

overload  cycle,  Aa  is  the  crack-size  increment  between  the  overload  and  current 
cycles,  and  z is  a parameter  which  defines  the  size  of  the  hardened  zone  ahead 
of  the  crack  tip.  The  latter  parameter  can  be  calculated  from  plastic-zone-size 
estimates : 


z (K  /f . ) 

^OL  max  ty 


(20) 


where  f is  the  material  tensile  yield  strength.  The  Willenborg  model  then 
ty 

defines  the  effective  stress  range  and  ratio  according  to: 


AK*  = K*  - K*  . = K “ 

max  min  max  R ^ 


- K ) = AK 
min  R 


(21) 


R*  = K*.  /K*  = (K  . - K )/(K  - K ) 

min  max  min  R max  R 


(22) 


Thus,  the  procedure  for  incorporating  Willenborg  retardation  in  the  load 
model  must  preserve  As  while  adjusting  the  stress  ratio. 


R = K . /K  = (S  - AS)/S 

min  max  max  max 


(23) 


If  it  is  assumed  that  K is  sufficiently  insensitive  to  local  geometry  effects, 
the  changes  in  these  effects  over  a few  load  cycles  can  be  neglected,  and 
Eq.  22  can  be  reformulated  in  terms  of  stress  parameters: 
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R* 


(24) 


= (S  - As  - 
max 


S^)/(S 
R max 


S^) 


where 


OL 


S = S /I  - Aa/z^ 

R max  OL 


(25) 


Once  s has  been  calculated  for  the  current  load  cycle,  R'  can  be  found  and 
R 

s can  be  adjusted  according  to: 
max 


5?  * = o _ 

max  1-R*  ^max  “ l^^T 


(26) 


The  forward  extent  of  the  retardation  effect  is  bounded  by 


Aa  < Aa  = /I  ^ , OL  2 

— max  OL  1 - (S  /S  ) 

max  max 


(27) 


Thus,  incorporation  of  Willenborg  retardation  in  the  load  model  is 
straightforward,  provided  that  a reasonable  method  for  estimating  Aa  can  be 
formulated.  One  possible  approach  is  to  include  cycle-by-cycle  da/dN  calcula- 
tions, using  the  Forman  equation,  in  the  cycle-by-cycle  inspection  of  loads. 

Since  these  calculations  extend  only  over  one  typical  set  of  (say)  50  flights 
per  mission  type,  the  added  computational  cost  is  not  excessive.  The  da/dN 
calculations  must  be  performed  in  parallel  with  several  initial-crack  sizes 
which  represent  the  entire  range  of  crack  size  which  might  be  expected  in  one 
airplane  life  (e.g.  MIL-A-83444  criteria  size  to  critical  size) . Obviously, 
the  retardation  effect  will  become  less  significant  late  in  an  airplane  life, 
when  the  current  crack  sizes  are  larger.  The  gradual  disappearance  of  retarda- 
tion will  be  reflected  in  the  statistics  of  the  adjusted  load  spectra  (i.e.  s*  , 

max 

as  given  by  Eg.  26)  , and  may  be  thought  of  as  a fictitious  usage  change  for  the 

purpose  of  random-load  crack-growth  prediction. 

Three  approaches  are  possible  for  the  loads-ad justment  procedure.  First, 

one  might  assume  that  the  effect  of  each  overload  extends  only  to  the  occurrence 

of  the  next  overload.  In  this  case,  the  adjustments  of  s fall  into  non- 
max 

overlapping  compartments.  Second,  one  might  assume  that  the  effect  of  each 
overload,  extends  to  the  occurrence  of  the  next  larger  overload,  in  which  case 
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the  compartments  overlap  partially.  Finally,  one  might  assume  that  retardation 
is  a cumulative  effect  which  applies  to  all  cycles  forward  of  a given  overload, 
subject  only  to  the  crack-growth  limitation  given  by  Eq.  27.  In  all  cases,  the 
adjustments  can  be  made  by  appropriately  controlled  scans  of  the  raw  reduced 
flight  data  As,  s . Selection  of  an  approach  and  determination  of  the  number 
of  representative  crack  sizes  required  to  permit  smooth  transitions  of  load 
statistics  over  an  airplane  life  are  subjects  which  require  additional  study. 

In  particular,  the  latter  determination  should  be  based  upon  the  results  of 
calculations  with  actual  flight-loads  spectra,  especially  since  the  amount  of 
detail  required  may  be  sensitive  to  usage. 

2. 3 Derivation  of  Crack-Size  Statistics 

The  analytic  prediction  of  crack  growth  statistics  was  considered  as  an 
alternative  to  analog  Monte  Carlo  simulation.  If  feasible,  it  would  eliminate 
the  necessity  of  utilizing  analog  computer  facilities.  At  the  very  least,  this 
analysis  would  result  in  a better  understanding  of  the  processes  involved,  so 
that  analog  simulation  could  be  performed  both  more  accurately  and  more  econom- 
ically. 

Analytic  prediction  of  the  statistics  proved  to  be  quite  easy  for  simple 
models  of  loading  and  crack  growth  behavior.  A FORTRAN  program  was  written 
to  predict  the  statistics  of  final  crack  size  for  a certain  type  of  crack- 
growth  problem,  specifically  the  problem  considered  in  previous  work  [19]. 

The  Paris  equation  is  used  to  model  crack  growth  behavior  and  the  loading  is 
modeled  as  follows.  The  total  number  of  loading  cycles,  N,  is  given,  and 
these  cycles  are  divided  evenly  into  n blocks.  Within  each  block,  the  value 
of  As  is  constant.  The  value  of  As  for  each  block  is  found  by  sampling  an 
input  noise  consisting  of  white  noise  and  an  additive  constant.  This  is  an 
extremely  crude  model  of  the  statistics  of  As,  expecially  since  it  even  pro- 
vides nonzero  probability  of  negative  As,  but  it  suffices  for  an  example. 

Since  the  Paris  crack-growth  model  was  utilized  with  = 4,  the  damage 
parameter  h was  chosen  as  h = 1/a,  with  the  rate  equation  as  given  by  Eg.  4. 
Since  the  rate  of  change  of  the  damage  parameter  does  not  depend  on  its 
instantaneous  value,  statistical  analysis  proceeds  with  little  difficulty. 

The  values  of  As  for  the  various  blocks  are  statistically  independent,  and 
they  have  the  same  given  normal  probability  density  function.  All  moments 
of  this  density  function  are  known,  where  these  moments  are: 
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(28) 


.(AS) 


% ~ E[(AS)  ] = Expected  value  of  (AS)^ 


Therefore,  the  moments  of  (As)  are  also  known,  being  given  by: 


’ = E[{  (AS)‘*}P]  = 

y 4p 


(29) 


The  total  damage  sustained  by  the  structure  after  the  N loading  cycles, 
with  constant  As  (constant  amplitude  cycles)  within  each  of  the  n blocks,  is  then 


(30) 


It  is  convenient  to  define: 


g = 


Ah 

NC 

p 


1 

n 


n 


(AS . ) '* 

1 


(31) 


which  will  eliminate  confusion  in  the  following  analysis.  The  expected  value 
of  g is  simply: 


E[g]  = g = E[(AS)  ] = 


4 

(AS^) 


(32) 


The  statistics  of  g are  represented  by  this  mean  value  and  by  the  central 
moments,  which  are  the  moments  of  the  probability  density  function  about  g. 
The  central  moments  are: 
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= E[(g-i)’^]  = 


n n 


n k 

. I (E[.n  {(As^  - g}]) 


(33) 


k 

In  general,  there  are  n terms  to  consider.  This  is  excessive,  and  is  reduced 
by  grouping  terms  of  the  same  type,  as  described  below.  The  total  effect  of 
each  group  is  calculated  and  added  to  the  effects  of  other  groups,  greatly 
reducing  the  number  of  computations  necessary. 

A "type"  is  distinguished  from  another  "type"  if  the  expected  value  of  the 
product: 


k 

E[  n 

j=l 


{(As.  )^  - g}] 
j 


(34) 


is  not  the  same.  A more  useful  method  of  distinguishing  between  types  can 
be  developed.  If  the  indices  i^  through  i^^  are  separated  into  groups,  with 
identical  values  grouped  together,  then  there  will  be  y groups  with  elements 
in  the  i^  group.  The  above  expected  value  can  then  be  written  as: 


k 

"fjSi 


{ (As 


,4  _ Y (As'^) 

g}]  kSi  ^(Y^) 


(35) 
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Therefore,  two  terms  are  of  the  same  type  if  they  contain  the  same  number  of 
groups  y,  and  if  the  values  of  are  the  same.  The  order  in  which  the  values 
of  appear  is  not  important. 

A procedure  has  been  developed  for  finding  all  possible  types  of  terms 
and  calculating  the  number  of  the  n terms  which  are  of  each  type,  without 
having  to  consider  the  terms  individually.  The  contributions  of  all  types 
are  combined  to  economically  provide  the  results  for  the  moments  of  g,  which 
can  easily  be  processed  to  provide  the  moments  for  Ah. 

As  the  number  of  blocks  (n)  increases,  the  probability  density  of  Ah  will 
approach  a normal  density  function.  All  odd  central  moments  will  approach  zero, 
and  all  even  central  moments  will  approach  the  appropriate  values  for  a normal 
density  function.  The  amount  by  which  a density  function  departs  from  the 
normal  can  be  represented  by  its  eccentricities,  which  are  defined  as; 


e„  = y /(y-,) 


n/2 


for  n odd 


(36) 


e = 
n 


n/2 


- Ix3x5x  ...  X (n-1) 


for  n even 


A probability  density  function  is  then  completely  described  by  its  mean, 
standard  deviation,  and  its  eccentricities  of  order  3 cind  higher.  Only  the 
first  few  eccentricities  are  necessary  for  a reasonable  approximate  representa- 
tion of  the  statistics,  since  only  the  first  few  are  likely  to  be  non-negligible. 

A FORTRAN  program  has  been  written  which  provides  the  mean,  standard  devia- 
tion, and  the  third  and  fourth  eccentricities  (skewness  and  excess  kurtosis)  of 
the  probcdDility  density  function  for  Ah.  This  program,  using  the  procedure  out- 
lined above,  is  presented  as  Appendix  A of  this  report.  Some  example  output 
from  the  program  is  presented  as  Appendix  B. 

2.4  Reanalysis  of  Previous  Monte  Carlo  Simulation 

A reanalysis  of  the  previous  work  [19]  on  analog  Monte  Carlo  simulation  of 
crack  growth  has  been  performed.  This  reanalysis  has  revealed  some  errors  in 
the  work,  and  hence  has  affected  the  details  of  the  continued  development  of 
the  Monte  Carlo  approach.  This  reanalysis  has  also  helped  to  compare  the  analog 
Monte  Carlo  approach  with  other  approaches  to  the  overall  problem. 
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The  previous  work  performed  Monte  Carlo  experiments  of  fatigue  crack 
growth  as  follows.  White  noise,  actually  limited  in  frequency  but  apparently 
containing  high  enough  frequencies  for  the  purpose,  was  produced.  The  noise 
with  an  added  constant  was  sampled  at  intervals,  and  the  value  of  each  sample 
was  held  until  the  next  sample  time,  resulting  in  a function  with  steps.  This 
function  was  assumed  to  represent  the  loading  parcimeter  As  in  the  Paris  equation. 
The  dynamics  of  crack  growth  as  modeled  by  the  Paris  equation  were  then  simulated 
on  the  analog  computer,  providing  crack  size  as  a function  of  time.  This  proce- 
dure was  rerun  many  times  for  the  same  initial  conditions,  with  varying  crack 
size  results  because  of  the  randomness  of  the  simulated  loading  parameter  As. 

The  varying  results  were  then  inspected  to  estimate  the  statistics  of  the  final 
crack  size  under  the  simulated  conditions. 

The  previous  work  used  the  Paris  crack  growth  model  with  = 4.  As  has 
been  described,  the  appropriate  choice  for  the  damage  parameter  in  this  case  is 
h = 1/a,  giving 


dN 


-C  (AS)"^ 
p 


(37) 


Consequently,  as  explained  in  Subsection  2.3,  experiments  with  the  same  loading 
but  with  different  initial  values  of  crack  size  a^  are  identical,  when  considered 
in  terms  of  the  change  in  the  damage  parameter  from  its  initial  value.  Ah  = h - h^. 

The  previous  work  consisted  of  two  classes  of  experiments.  Both  had  As 
modeled  as  sampled  white  noise  with  standard  deviation  and  an  added  constant 

Xs.  Within  each  class  these  parameters  were  the  same,  but  initial  crack  size  was 
different  for  the  different  experiments.  Each  class  can  then  be  considered  as  a 
number  of  identical  experiments  in  Ah,  with  different  scaling  factors  used  in 
modeling  the  experiments  on  the  analog  equipment. 

The  experimental  data  to  be  reanalyzed  were  available  on  punched  cards.  The 
data  were  re-interpreted  in  terms  of  Ah  providing  the  following  four  results  for 
each  experiment: 

mean  of  Ah 

standard  deviation  of  Ah 
3^  = skewness  of  Ah 
(3^  - 3)  = excess  kurtosis  of  Ah 

22 


Ah  = 


^Ah 


fAh) 


(Ah) 


The  skewness  and  excess  kurtosis  represent  the  amount  by  which  the  probability 
density  function  for  Ah  differs  from  a normal  distribution.  The  above  four 
parameters  were  computed  for  four  times  in  the  life  of  each  experiment  of  the 
previous  work  (8K,  lOK,  13. 3K,  and  20K  cycles). 

The  means  Ah  were  found  to  be  fairly  consistent  within  each  class  of 
experiments.  The  standard  deviations  were  less  consistent,  with  variations  of 
20  percent  between  high  and  low  values  not  uncommon.  The  third  and  fourth 
eccentricities  (skewness  and  excess  kurtosis)  were  sometimes  more  thain  an  order 
of  magnitude  apart.  Considering  that  each  experiment  provided  about  200  reruns 
for  the  data  analysis,  the  experimental  procedure  must  be  seriously  questioned. 

The  previous  work  reported  that  as  the  sampling  rate  was  varied  up  to 

about  200  Hz,  the  standard  deviation  of  Ah  was  reduced.  However,  as  the  sampling 

rate  was  increased  above  200,  even  to  1200  Hz,  the  standard  deviation  did  not 
continue  to  decrease.  This  is  in  conflict  with  the  known  results  described  in 
the  previous  subsection,  and  is  therefore  another  reason  to  seriously  question 
the  experimental  procedure.  Because  of  this  abnormality,  we  attempted  to 
reproduce  the  previous  work,  utilizing  the  same  electronic  equipment.  The 
sample-and-hold  circuitry  was  immediately  recognized  as  the  primary  cause  of 
the  discrepancies. 

Figure  7 shows  some  examples  of  sampling  of  white  noise  using  the  solid- 
state  sample-and-hold  mechanism  on  the  analog  computer.  Although  the  white- 
noise  generator  appears  to  be  adequate,  the  samples  are  obviously  correlated 

(each  sample  is  likely  to  be  similar  to  the  immediately  preceding  sample) . It 

is  known  that  samples  of  white  noise  are  completely  uncorrelated  regardless  of 
the  sampling  rate.  This  explains  the  observed  discrepancies  in  the  experimental 
results. 

The  prediction  method  for  crack  growth  statistics  described  in  Section  2.3 
was  applied  to  the  two  classes  of  experiments.  All  of  the  input  parameters  for 
the  calculation  were  known  except  for  the  number  of  samples.  Sampling  was  performed 
at  250  Hz,  which  would  provide  (at  the  four  times  considered)  40,  50,  66.5,  and 

100  samples.  However,  because  of  the  correlations  of  the  samples  introduced  by 

the  sample-and-hold  mechanism,  the  effective  sampling  rate  was  much  slower.  If 

the  calculations  are  performed  for  approximately  one  fifth  the  above  sample  rate 

(9,  11,  15,  and  22  samples  at  the  four  times  considered),  then  the  experimental 
and  analytical  results  are  fairly  consistent.  The  mean  Ah  is  consistently  about 
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10  percent  higher  experimentally  than  analytically,  and  since  it  does  not 
depend  on  sampling  rate  this  discrepancy  is  attributed  to  other  analog  circuit 
errors.  Comparison  of  values  of  also  shows  good  agreement,  but  of  course 
this  agreement  is  caused  by  the  choice  of  sampling  rate  used  in  the  analysis, 
one- fifth  the  actual  rate.  The  conclusion  to  be  drawn  from  this  is  that  the 
sample- and- ho Id  mechanism  introduced  a correlation  time  of  about  five  samples 
into  its  output.  Inspection  of  Fig.  7 shows  that  this  correlation  time  is 
reasonable. 

To  summarize,  it  must  be  concluded  that  the  previously  reported  statistics 

of  final  crack  size  were  primarily  the  result  of  a spurious  effect  in  the  analog 

sample-and-hold  device,  which  was  therefore  rejected  for  subsequent  work.  The 

problem  with  the  sample-and-hold  device  was  unfortunate,  but  not  serious,  since 

the  device  had  only  been  utilized  as  a temporary  measure  pending  more  accurate 

modeling  of  the  loading  characteristics.  The  more  accurate  modeling  has  been 

accomplished  by  the  interim  model  for  time-histories  of  s . and  s mentioned 

min  max 

in  Subsection  2.2. 

2.5  Utilization  of  Density  Moments 

As  previously  described,  the  analytic  derivation  of  final  crack  size 
statistics  results  in  the  first  few  moments  of  the  probability  density  function 
of  a damage  parameter,  h.  Once  these  moments  have  been  determined,  they  must 
be  used  to  develop  a model  of  the  probability  density  function  itself.  This 
will  then  allow  calculation  of  probabilities  of  interest,  such  as  structural 
failure  or  crack  size  large  enough  for  visual  detection.  These  probabilities 
of  interest  are  mathematically  represented  by  inequalities  in  the  damage 
parameter. 

It  is  not  difficult  to  formulate  models  of  probability  density  functions 
which  have  the  specified  first  few  moments.  It  is  much  more  difficult  to  formu- 
late such  models  which  are  realistic.  The  first  attempt  was  formulated  as; 

2 , n-i 

f(h)  = f(h)  {1  + + a^h  + a^h  + ...  + a^h  ) os) 

where  f (h)  is  the  density  function  for  h,  and  f(h)  is  a normal  density  function 
with  the  same  mean  and  standard  deviation.  Unfortunately,  although  this  exactly 
matches  the  specified  first  n moments  of  the  density  function,  the  unspecified 
higher  moments  were  extremely  unusual.  The  resulting  model  of  the  density 
function  appears  valid  near  the  mean,  but  far  from  the  mean  the  function  has 
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negative  values  which  are  physically  unrealizable.  This  approach  would  there- 
fore be  sufficient  for  computation  of  exceedance  probabilities  near  the  mean. 
However,  more  extreme  values  cannot  be  evaluated  by  this  approach. 

Consequently,  another  approach  must  be  taken.  The  second  approach 
considers  the  log-characteristic  function  of  the  probability  density, 

=1X1  [ r f(h)dh] 

(39) 


where  j = This  is  simply  the  logarithm  of  the  inverse  Fourier  transform 

of  the  probability  density  function.  For  a normal  density,  this  function  is 
simply  a quadratic. 


1(^(0)) 


- 1 2,  2 
= jb)h  - — U)  h 


(40) 


The  model  to  be  used  for  crack-size  statistics  incorporates  a higher-order 
polynomial  in  h,  which  will  satisfy  specified  values  of  the  first  few  moments, 
while  giving  reasonable  values  for  the  unspecified  higher  moments  of  the 
probability  density  of  h.  This  approach  has  been  used  previously  for  similar 
problems  [27]. 

For  example,  when  the  first  four  moments  of  a random  variable  are  specified 
as 


m 


1 


h 


(41) 


and  the  higher-order  moments  are  unspecified,  the  log-characteristic  function 
ip  (o))  would  be  a quartic  in  O),  and  the  first  four  moments  of  h are  exactly 
matched  if  [28]: 


j + 


1.3,  2 

- — DOJ  (m^  - + 2mj^)  + 


+ 3^  oj'*  (m^  - 4m3m^  - + L2m^  - emj) 
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(42) 


However,  the  computation  of  the  log- characteristic  function  is  only  a 
first  step.  The  desired  results  are  exceedance  probabilities  for  extreme 
values  of  h,  which  will  be  important  in  risk  analyses.  This  is  theoretically 
accomplished  by  applying  a Fourier  transformation  to  the  exponential  of  \\)  (co) 
to  provide  the  probability  density  f (h)  and  then  taking  appropriate  semi- 
infinite integrals  of  f (h)  to  provide  the  exceedance  probabilities.  Neither 
step  can  be  accomplished  analytically;  numerical  methods  will  be  required. 
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Section  III 


ANALOG  SIMULATION  OF  CRACK  GROWTH 


3.1  Results  for  Specific  Tasks 

Eight  of  nine  specific  analog  programming  tasks  were  completed  in  the 
present  investigation.  The  following  subsections  outline  each  task,  describe 
the  procedures  used  to  implement  the  task  on  the  analog  computer,  and  present 
the  results  of  test  runs.  For  the  purpose  of  the  study,  the  simulations  of 
all  expressions  for  stress  intensity  factors  have  been  simplified  by  omitting 
the  constant  factors,  i.e.  K = s/rra  for  a crack  in  an  infinite  plate  is 
replaced  by  K ^ s/a.  (in  any  case,  the  constant  factors  would  have  been 
absorbed  in  the  Forman  equation  rate  constant,  All  test  runs  were  made 

on  the  Cl- 5000/5  hybrid  analog/digital  computer  at  the  USAF  Aeronautical  Systems 
Division  (ASD/ADSD) , Wright-Patterson  Air  Force  Base,  Ohio. 

3.1.1  Forman  Equation 

Simulation  of  the  Forman  equation  on  the  analog  computer  is  a straight- 
forward task.  The  simulation  was  programmed  in  terms  of  the  inverse  damage 
parameter  "h"  (see  Eqs.  12-14  in  Subsection  2.1).  For  the  purpose  of  this 
study,  7075-T6  aluminum  was  taken  as  an  example  alloy,  with  experimental 
growth-rate  data  [20]  fit  by: 

= 3 C^  = 5x10”^^  = 68  ksi/in. 


For  the  special  case  of  N = 3,  Eq.  13  can  be  simplified  to  the  form: 

F 


dN 


-C  (AS)  ^ 


S 

max 


S 

c 


- S 

max 


(44) 


where  S = K h and  C = (1-N  /2)C  /K  . 

G c F F c 

The  analog  flow  chart  for  Eq.  44  appears  in  Fig.  8.  Stability  and 

accuracy  were  checked  initially  by  a simulation  with  constant- amplitude 

loading  (s  . and  s held  constant,  cycles  converted  to  analog  computer 
min  max  2 

time  by  analog  scaling  laws).  An  initial  condition  of  h^  = 10  (a^  = 1/h  = 0.01 

inch)  was  specified.  The  solution  was  stable  until  a crack  size  of  approximately 
(h  = l//a  = 0.5)  was  reached. 
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Typical  results  of  additional  tests  with  random  load  inputs  are  illustrated 

in  Figs.  9 and  10.  Band-limited  white  noise  with  a 350  Hz  cutoff  frequency  was 

used  to  model  both  s . and  s . These  signals  are  fully  correlated  because 

min  m^x 

they  were  produced  from  the  same  noise  generator  without  time-delay.  Hence, 

As  is  still  a constant-amplitude  signal,  as  is  evident  in  the  figures.  However, 
dh/dt  is  still  a random  variable,  since  s^^^  appears  alone  in  Eq.  44.  No  over- 
loads occurred  in  the  analog  operational  amplifiers  when  the  unfiltered  noise 
input  was  used.  Figure  9 presents  a single,  real-time  solution.  Figure  10 
illustrates  a sequence  of  solutions  with  the  analog  run  in  its  repetitive 
operation  mode.  The  several  replicates  run  in  this  mode  resulted  in  identical 
values  of  final  crack  size,  the  value  being  in  agreement  with  analytical  calcula- 
tion based  on: 


(dh/dN)  = -C  (AS)^  S 

lUaX 


c max 


(45) 


This  is  in  accordance  with  the  property  of  /dh/dN  (based  on  the  Central  Limit 
Theorem)  discussed  in  Section  2. 

3.1.2  Retardation 

Retardation  models  such  as  the  Wheeler  model  [14]  do  not  account  directly 
for  load- interaction  effects.  The  Wheeler  model  might  be  applied  to  the  Forman 
equation  by  replacing  Eq^  9 with: 


da 

dN 


^F  , /f  ^f/2  m ^max^  . 
— [(AS)  a ] [-^Hg 

c 


- s 

c max 


(46) 


The  Wheeler  exponent,  m,  can  be  determined  by  trial-and-error  adjustment  until 
acceptable  predictions  are  obtained  for  crack  size  versus  cycles  by  comparison 
with  experimental  data  for  a given  load  spectrum.  It  is  evident  that  Eq.  46 
can  be  simulated  with  no  more  difficulty  than  encountered  in  simulating  a gen- 
eral Forman  equation.  One  need  only  define  an  appropriate  damage  parameter  and 

define  S accordingly.  For  instance,  for  mN  > 2,  the  "natural"  definition  of 
c ^ 

the  damage  parameter  is : l-mNp/2 

h = a (47) 
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cind  S is  then  found  to  be 
c 

S = K /v^  = K (48) 

c c c 

Hence,  the  only  extra  step  for  a Wheeler-Forman  simulation  involves  the 
determination  of  the  exponent,  m,  from  laboratory  tests. 

However,  the  currently  preferred  interaction  models  account  individually 
for  the  effects  of  isolated  overloads.  The  present  program  task  required  a 
demonstration  of  the  ability  to  simulate  retardation  by  isolated  overloads, 
including  a gradual  return  to  unretarded  growth  rates.  If  one  of  the  accepted 
interaction  models  [15,16]  is  adopted,  the  simulation  then  implies  that  point- 
calculations  of  the  retardation  effect  must  be  made  from  continuous- function 
representations  of  the  stress  parameters.  Such  a procedure  would  undoubtedly 
increase  the  analog  computing  error  to  unacceptable  levels. 

Additional  insight  into  the  complications  involved  can  be  gained  by 
considering  a Paris  equation,  modified  by  a simplified  interaction  model. 
Suppose  that 

N 

Pq 

da/dN  = Cp(As/a)  (49) 

in  the  absence  of  overloads,  but  that; 

N* 

da/dN  = C (AS»/a)  ^ (^0) 

P 


N' 

P 


N (1-e 


■An/n 


) + 


-An/n 
N e 

Pi 


(51) 


replaces  Eq.  49  immediately  after  an  overload.  In  Eqs.  50  and  51,  the  cycle 
count  An  begins  at  zero  on  the  cycle  immediately  following  the  overload,  and 
N , N are  assumed  to  be  given  quantities  with  N < N . Hence,  the  retarda- 

Pi  y 

tion  effect  appears  as  a series  of  step-drops  followed  by  exponential  rises  in 
the  growth- rate  exponent.  One  can  conceive  of  simulating  this  model  on  the 
analog.  However,  there  is  no  longer  any  unique  relation  between  crack  size 
and  the  damage  parameter: 


h 


1-Ny2 
a P 


(52) 
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Therefore,  a simulation  directly  in  terms  of  crack  size  would  be  required. 
Furthermore,  the  principal  advantage  of  the  analog  (using  a short  period  of 
computing  time  to  represent  a large  number  of  load  cycles)  is  lost  because  of 
the  need  to  represent  the  individual  overload  cycles.  Thus,  if  time-compression 
is  used,  many  isolated  overloads  appear  lumped  together  in  the  stress-parameter 
signals,  and  Eq.  50  is  activated  too  infrequently.  Similar  conclusions  can  be 
drawn  about  simulation  of  the  Forman  equation  with  an  accepted  retardation 
model.  Therefore,  retardation  models  were  not  implemented  on  the  analog  computer. 

3.1.3  Local  Geometry  Effects 

The  third  programming  task  was  to  demonstrate  the  ability  to  vary  a 
simulated  crack  growth  rate  by  including  a local  geometry  effect,  e.g.  the 
change  in  stress  intensity  factor  for  a crack  which  grows  away  from  a fastener 
hole.  The  geometry  effect  can  be  expressed  in  terms  of  a dimensionless 
parameter,  3,  such  that  the  effective  range  in  stress  intensity  factor  is 
given  by; 


AK 


eff 


BAk 


(53) 


The  parameter,  3,  can  usually  be  stated  in  terms  of  the  ratio  of  crack 
size  to  a significant  dimension  associated  with  the  geometry  of  a local 
structural  detail.  However,  3 cannot  always  be  expressed  analytically.  For 
example,  3 is  a function  of  a/r  for  a crack  growing  away  from  a fastener  hole, 
where  r is  the  radius  of  the  hole.  If  the  hole  is  open  and  is  located  in  a 
large  panel  subjected  to  uniform  edge  tension,  and  if  the  crack  orientation 
is  perpendicular  to  the  tension,  then  3=  3.4  for  very  small  cracks  (a/r  0) 
and  3 I//2  as  a/r  «>,  according  to  a complex  variable  solution  by  Bowie  [29]. 

This  and  many  other  numerical  solutions  for  3 associated  with  other  simple 
geometries  are  available  in  graphical  form  [30] . For  more  complicated 
geometries,  3 can  be  obtained  by  finite-element  analysis,  e.g.  using  hybrid- 
stress  crack-tip  elements  [31]. 

If  the  local  geometry  effect  is  included  in  the  Forman  equation,  the 
crack  growth-rate  becomes; 


da 

dN 


""f 

c^(6Ak) 

r 

(l-R)K  - BAK 
c 


(54) 
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Note  that  R = K . /K  is  not  affected  by  6.  Equation  54  can  then  be 
min  max  ^ 

transformed  into: 


dN 


» - ^ r -s  -'es  ■ 

c c max 


(55) 


where  Ak  ~ As/a  has  been  assumed  as  before.  In  principle,  6 can  be  described 
directly  in  terms  of  the  damage  parameter,  rather  than  crack  size. 

In  the  present  study,  the  empirical  finite-width  correction  factor: 


3=1+  1.1896X^  + 1.3016X^  + 1.3650X®  + 
+ 1.3739A®  + 1.4764A^° 


(56) 


was  chosen  for  implementation.  Equation  56  represents  the  accelerating  effect 
of  a single  edge-crack  of  length,  a,  growing  across  a plate  of  width,  b,  with 
the  plate  subjected  to  uniform  tension,  and  with  X = a/b.  Thus,  one  could 
calculate  X continuously  on  the  analog  by  substituting  the  damage-parameter 
transformation  to  obtain: 


^ ^ ^2/(2-Np) 


(57) 


However,  a function-generation  approach  is  obviously  more  convenient 

and  more  economical  in  terms  of  the  number  of  analog  components  used.  Most 

modern  analog  computers  include  function  generators  which  nominally  provide 

piecewise  linear  representations  for  functions  of  one  input  variable,  with  ten 

to  twenty  segments  allowed.  Equation  56  provides  a severe  test  of  this  capability 

because  of  the  strong  nonlinearity  in  3.  Figure  11  compares  Eq.  56  with  a typical 

representation  obtained  on  the  analog.  In  this  case,  the  simulation  has  been 

formulated  in  terms  of  1/3  to  avoid  numerical  instability  in  the  function  generator. 

The  points  were  selected  for  segment  generation  for  a plate  width  b = 10  inches 

and  rate  exponent  N = 3,  as  follows: 

F 
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Crack  Size,  a 

(inch) 

Function  Input, 
h=l/v^  (in.  ^^^) 

Function  Output, 

1/6 

0.01 

10 

1.0000 

0.01562 

8 

1.0000 

0.04 

5 

1.0000 

1.0 

1 

0.9743 

1.235 

0.9 

0.8797 

1.562 

0.8 

0.5709 

2.041 

0.7 

0.2010 

2.778 

0.6 

0.0467 

4.0 

0.5 

0.0078 

10.0 

0.316 

0.0000 

The  above  table  suggests  that  the  crack-growth  simulation  would  probably  have  to  be 
restricted  to  a maximum  crack  size  between  1 and  2 inches  for  this  case.  Figure  11 
illustrates  the  limitation  of  the  analog's  ability  to  match  a severely  nonlinear 
function. 

With  regard  to  implementation  of  a geometry  factor  in  the  Forman  equation, 

it  is  apparent  that  B can  be  inserted  in  the  appropriate  places  in  the  flow  chart 

shown  in  Fig.  8,  with  the  aid  of  several  additional  multipliers.  One  multiplier 

is  required  to  form  the  product  3 s just  before  the  inverter  which  outputs 
^ ind.x  ^ 

S - s while  three  multipliers  are  required  to  form  3 and  then  the  product 

c max 

3^  As^  near  the  upper  center  of  the  diaqram.  The  solution- instability  problem 
will  not  exist  for  other  types  of  geometry  factors,  such  as  3 for  a crack  growing 
away  from  a fastener  hole.  Also,  the  exact  function  will  be  easier  to  match 
because  the  nonlinearity  is  much  less  severe. 

3.1.4  Randomized  GAG  Sequence 

For  the  purpose  of  predicting  crack  growth  a GAG  cycle  can  be  treated  as 
having  two  components,  both  of  which  result  from  the  sudden  change  in  mean  stress 
as  an  airplane  takes  off  or  lands.  For  example,  when  the  airplane  taxies  and 
during  the  pre- takeoff  roll,  a point  on  the  lower  skin  of  the  wing  experiences 
small— amplitude  stress  oscillations  about  a negative  (compression)  mean.  As 
the  airplane  approaches  takeoff  speed,  the  wings  bear  an  aerodynamic  lift  which 
approaches  the  aircraft  weight  and  reverses  the  sign  of  the  wing  static  bending 
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moment.  Following  takeoff,  the  wing  lower  skin  experiences  somewhat  larger 

stress  oscillations  about  a positive  mean.  One  can  thus  identify  two  effects 

on  crack  growth.  First,  the  initial  change  in  mean  stress  is  equivalent  to  a 

single,  relatively  large  range  at  the  beginning  of  each  flight.  The  magnitude 

of  the  GAG  range  can  be  treated  as  deterministic  (effect  of  mean  only) , as  a 

deterministic  range  from  last  taxi  minimum  to  first  gust  maximum,  or  as  random 

with  a small  variance  equal  to  the  sum  of  the  variances  associated  with  taxi 

and  gust  loads.  The  second  effect  is  that  the  change  in  mean  stress  must  be 

accounted  for  in  the  flight- loads  cycles.  This  effect  of  the  mean  change  on 

the  two-parameter  load  model  is  that  the  values  of  s increase  after  takeoff, 

max 

while  the  values  of  As  are  not  directly  affected  by  the  GAG  event. 

The  fourth  programming  task  required  a demonstration  of  the  ability  to 
randomize  the  sequence  and  time  of  appearance  of  a limited  number  of  GAG  ranges 
having  deterministic  magnitudes.  This  task  was  easily  accomplished  by  using 
band- limited  zero-mean  Gaussian  noise  as  a switching  input  summed  with  an 
adjustcible  bias  voltage,  -a,  to  make  the  signal  mean  negative.  With  an 
electronic  comparator  set  to  switch  "on”  when  the  total  switching  voltage 
exceeds  zero,  a constant  voltage  As^  can  be  added  intermittently  to  the 
continuous  signal  As(t)  which  represents  flight-loads  stress  ranges.  Figure  12 
illustrates  the  required  analog  flow  chart. 

The  average  number  of  inputs  of  As  per  second  can  be  estimated  from  the 

o 

expected  number  of  upward  crossings  of  zero  by  the  switching  signal  [32]: 


e(0) 


2 2 

1 - -a  /2a 
— V e 

2 O 


(58) 


where 


0^  = ^“’G(f)df  (59) 

and  where  G(f)  is  the  power  spectral  density  of  the  Gaussian  noise  in  units  of 
2 

(volts)  /Hz.  For  band-limited  white  noise,  Eqs.  59  and  60  are  easily  evaluated, 
resulting  in: 
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( 


(61) 


V = 2f 
0 c 


where  is  the  cutoff  frequency.  The  value  of  a can  then  be  determined  from 

Eq.  58  to  provide  the  proper  average  time  between  GAG  cycles. 

Figure  13  illustrates  some  typical  results  from  a series  of  test  runs  in 

which  As(t)  = 0 to  isolate  the  effect  of  As  . The  bias  voltage,  a,  was  varied 

o 

from  -1.51  to  -0.91  volts,  with  As^  = 10  volts.  The  Gaussian  noise  switching 
input  is  shown  as  x(t),  and  has  1 volt  rms  with  a 350-Hz  cutoff  (0  = 1, 

= 405  sec  ^)  . 

The  increase  of  average  time  between  GAG  cycles  with  decreasing  a is 

evident.  The  GAG  cycle  values  appear  to  be  random  only  because  the  pen-chart 

recorder  response  is  not  rapid  enough  to  reproduce  As^  faithfully.  The  actual 

As^  can  easily  be  varied  over  a small  range  to  represent  the  random  effect  of 

first-gust-peak/ taxi-minimum  by  adding  a small-amplitude  low-frequency  sinusoidal 

signal  to  As  . 

o 

3.1.5  Superposition  of  Random  Flight  Loads  and  GAG  Cycles 

The  fifth  programming  task  required  a demonstration  of  the  ability  to  simulate 
the  second  effect  of  the  GAG  cycle,  i.e.  the  change  in  mean  stress  (and  s ) 
following  takeoff  and  landing.  The  implementation  of  this  task  requires  two 
comparators  and  a time  clock.  The  time  clock  must  be  synchronized  with  the  start 
of  the  problem  solution.  Therefore,  one  simply  arranges  to  compute; 


t* 


(62) 


on  the  analog,  where  t represents  the  current  analog  time  and  C is  a scaling 

factor  chosen  for  convenience  (to  prevent  the  clock  from  saturating) . 

If  the  two  comparators  are  fed  with  constant  switching  voltages  t!^  and 

t^  > tj^,  they  can  be  used  to  choose  between  three  load  environments,  with  the 

selection  changed  at  computing  times  t = t^/C,  t^/C.  The  analog  circuit  shown 

in  Fig.  14  provides  such  a choice  between  three  random  signals  which  have  (in 

general)  different  means  and  different  variances.  The  output  signal  s (t)  might 

then  be  used  to  simulate  s . Figure  15  illustrates  the  performance  of  this 

max 

circuit  with  the  analog  in  repetitive  operation  mode., 
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3.1.6  Random  Initial  Condition 


The  sixth  prograunming  task  required  a demonstration  of  the  ability  to 
simulate  a suitable  distribution  of  initial-crack  sizes.  Such  a distribution 
would  be  sampled  infrequently  (once  per  simulated  airplane  life)  to  provide 
the  initial  conditions  for  a Monte  Carlo  simulation.  A non-Gaussian  model  is 
desired  because  initial-crack  size  distributions  are  usually  skewed. 

This  task  is  implemented  by  using  an  appropriately  shaped  zero-mean 
white-noise  source.  The  shaping  function  is  derived  simply  as: 

A(x)  = f(x)/4)(x)  (63) 


where 

X = x(t)  = Noise  signal 
A(x)  = Required  shaping  function 
f (x)  = Sought- for  probability  density  of  output 
signal  amplitudes 

(j)  (x)  = Probability  density  of  input  noise  cimplitudes 

In  the  present  case,  cf)  (x)  corresponds  to  a Gaussian  source  with  prescribed 
rms  voltage,  0, 


2 2 

/ N 1 -x^/2a^ 

a/27r 


— oo  < X < 


(64) 


while  a two-paraimeter  Weibull  distribution  was  sought  for  the  output  noise: 


0 < X < CO  ; > 0 

b“ 


(65) 


Hence  the  required  shaping  function  becomes: 

ow/irr 

3 20 


A(x)  = IxP  " exp 


(66) 


Note  that  the  shaping  function  must  be  symmetric  about  x = 0 to  account  for  the 
restriction  of  the  output  to  positive  values. 
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Figure  16  compares  the  required  shape  with  a gain  function  actually 
obtained  from  one  of  the  analog  function-generator  components.  The  required 
shape  is  intended  to  produce  a Weibull  distribution  with  a = 2.5  and  3 = 0.007 
from  Gaussian  noise  with  O = 0.1.  This  is  another  case  where  the  analog 
component  has  difficulty  in  following  a severely  nonlinear  function.  (The 
gain  function  shown  was  obtained  only  after  seven  function-generator  circuits 
had  been  tested  for  accuracy  and  rejected.) 

The  eighth  gain  function  obtained  (Fig.  16)  was  given  an  accuracy  test  by 
sampling  the  noise  signal  1,000  times  at  sampling  rates  of  50,  100  and  200  Hz. 
The  1,000  sample  values  in  each  case  were  then  used  to  obtain  maximum- likelihood 
estimates  for  the  output  distribution  parameters,  according  to  [33] : 
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where  N is  the  total  number  of  samples.  The  results  obtained  were  as  follows: 


SAMPLING 

MAXIMUM-LIKELIHCX)D 

ESTIMATES  FOR: 

RATE  (Hz) 

Shape  Parameter,  a 
(Required  Value=2.5) 

Scale  Parameter,  3 
(Required  Value=0.007) 

50 

2.01 

0.0073 

100 

2.63 

0.0071 

200 

2.84 

0.0069 

The  above  results  indicate  that  there  may  be  some  difficulty  in  maintaining 
the  distribution  shape  parameter  at  the  much  lower  sampling  rate  which  would 
correspond  to  one  simulated  airplane  life. 

3.1.7  Simulation  of  Random  Maneuver  Loads 

The  production  of  random  noise  to  simulate  stress  time-histories  is  a task 
similar  to  the  simulation  of  a random  population  of  initial-crack  sizes.  The 
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programming  task  requires  a simulation  of  a non-Gaussian  process,  such  as 
might  represent  fighter  maneuvering  loads. 

In  this  case,  it  was  decided  to  produce  Rayleigh  noise: 

2 2 

f(x)=^e^  ^ ; 0£x<oo  (69) 

Equation  69  also  describes  a Weibull  distribution  for  the  special  case  a = 2, 
g = a/2.  The  shaping  function  required  to  produce  Rayleigh  noise  from  zero- 
mean  Gaussian  noise  is: 


A(x) 


(70) 


Hence,  the  shaping  filter  can  be  simply  constructed  with  a constant  gain,  two 
inverters  and  a diode,  as  shown  in  Fig.  17. 

A sampling  experiment  similar  to  the  one  described  in  Subsection  3.1.6 
was  conducted,  using  a 1-rms  Gaussian  source,  with  the  following  results: 


SAMPLING 

MAXIMUM-LIKELIHOOD  ESTIMATE  FOR: 

RATE  (Hz) 

Shape  Parameter,  a 
(Required  Value=2) 

Scale  Parameter,  g 
(Required  Value=/2) 

50 

1.92 

1.43 

100 

2.04 

1.50 

200 

2.18 

1.54 

In  this  case,  the  shape  parameter  appears  to  be  less  sensitive  to  the  sampling 
rate,  probably  because  of  the  linearity  of  the  shaping  function. 

3.1.8  Incorporation  of  Threshold  Effect 

The  eighth  programming  task  required  incorporation  of  logic  circuitry  in 
the  Forman  equation  to  recognize  when  AK  is  less  than  a prescribed  threshold 
value  below  which  crack  growth  is  assumed  not  to  occur.  When  the  Forman 
equation  is  employed,  the  threshold  value  is  commonly  assumed  to  depend  upon 
the  stress  ratio,  e.g.: 

- (1-R)  <^1) 

where  AK^  is  a material  property  measured  in  a zero- to- tension  (R  = 0)  crack- 
TH 

propagation  test.  However,  the  stress-ratio  effect  on  threshold  was  ignored 
in  the  present  study,  i.e.  AK^^  = Ak^^  was  assumed. 


37 


The  threshold  effect  is  easily  implemented  by  placing  an  electronic 

comparator  between  the  gain,  C , and  the  integrator  which  outputs  the  damage 

parameter,  -h  (lower  left  in  Fig.  8).  With  the  switching  function  As/h  - 

supplied  to  the  comparator,  the  input  to  the  integrator  can  be  switched  between 

Eq.  13  when  the  threshold  is  exceeded  and  zero  when  AK  is  below  the  threshold. 

A partial  flow  chart  for  the  additional  circuitry  is  illustrated  in  Fig.  18. 

An  example  real-time  simulation  of  the  Forman  equation  with  the  threshold 

effect  was  run  with  loading  in  which  s . was  held  constant,  while  s was 

min  max 

modeled  by  white  noise  to  randomly  trigger  the  threshold.  The  results  are 

shown  in  Fig.  19,  which  records  the  switch  output  and  AK(t).  It  is  evident 

that  the  rate  equation  spends  more  time  in  the  circuit  as  crack  size  increases, 

until  virtually  all  AK  values  contribute  to  crack  growth.  Also,  it  is  apparent 

that  an  equivalent  simulation  can  be  implemented  with  Ak  as  given  by  Eq.  71, 

TH 

since  the  term  1-R  can  be  produced  as  s /As  using  one  additional  "X/Y” 

^ max 

component. 

3.1.9  Coupled  Dynamics  for  Two  Cracks 

The  final  programming  task  required  a demonstration  of  the  ability  to 
simulate  the  simultaneous  growth  of  two  or  more  cracks  with  the  rate  equations 
coupled  through  change  of  stress  intensity  near  one  crack  as  a function  of  the 
size  of  the  others.  Situations  of  this  type  may  occur  in  aircraft  structures, 
e.g.  for  two  cracks  180®  apart  growing  away  from  the  same  fastener  hole. 

In  the  present  study,  a simplified  situation  was  formulated  strictly  for 
the  purpose  of  demonstrating  the  coupling  effect  without  introducing  the  compli- 
cations associated  with  local  geometry  factors.  Suppose  that  two  panels  in 
parallel,  of  widths  and  L^,  are  loaded  such  that  both  panels  are  subjected 
to  the  same  elongation.  Let  the  panels  contain  edge-cracks  of  lengths  a^  and 
a^f  and  assume  for  the  purpose  of  the  demonstration  that  the  panel  compliances 
are  concentrated  in  narrow  strips  near  the  cracks.  Then  the  stress  levels  in 
the  two  panels  are  given  approximately  by: 


S 

1 


P(L^-ai)/Li 


(L^-a^)ti  + {Va2)t2 


(72) 


®2  (L-a,  )t,  + (L-a_)t, 


1 1 


2 2 


(73) 


38 


where  P is  the  total  load  applied  to  the  two  panels  and  where  t^,  t^  are  the 
panel  thicknesses.  The  problem  is  obviously  artificial,  but  still  serves  to 
test  the  ability  of  the  analog  to  simulate  coupled  Forman  equations. 

For  the  case  = L,  ~ ^2  ~ ^ ~ rate  equations  can 
be  expressed  as: 


dh. 
1 

dN 


"i>'i 

c 


- P F. 
c . max  1 
1 


];  i=lr2 


(74) 


where 


F.  (h  . 
1 1 


(L  - l/hJ)/Lt 

' 2 2 ' 

2L  - lA^  - IA2 


i=l,2 


(75) 


The  loading  is  described  by  the  parameters  P and  AP  = P - P . . The 

analog  flow  chart  appears  in  Fig.  20. 

Debugging  of  this  simulation  was  difficult,  due  to  the  instabilities 
2 

caused  by  the  1/h.  terms  in  the  equations.  However,  two  simulations  were 

^ -13 

successfully  completed  using  7075-T6  properties  (N  = 3,  C = 5 x 10  , 

r r 

K = 68  ksi  /in. ) . In  the  first  simulation  (Fig.  21) , the  initial  conditions 
c 

were  0.01-  and  0.01888-inch  cracks  (h^  = 10,  h^  = 7.3).  (The  values  of  the 

constants  C' , C'  are  equal  to  C /2K  times  the  analog  scale  factor.)  In  the 
12  F c 

second  simulation  (Fig.  22),  the  initial  conditions  were  equal-sized  cracks, 

while  the  growth  rate  of  the  second  crack  was  arbitrarily  increased  by  a 

-12 

factor  of  three  (C  = 1.5  x 10  ).  In  both  cases,  the  panel  dimensions 

^2 

= 10  inches  and  t^  = t^  = 0.1  inch  were  used,  while  constant- amplitude 

loading  with  AP  = P = 140  x 10^  lbs.  was  chosen.  In  each  case,  the  crack 
max 

which  grows  faster  is  observed  to  decelerate  toward  the  end  of  life,  while  the 
slower  crack  accelerates.  This  effect  results  from  the  transfer  of  load  from 
the  panel  with  the  larger  crack  to  the  panel  with  the  smaller  crack  because  of 
the  difference  in  panel  compliance. 

3.2  Evaluation  of  Results 

The  results  of  the  simulation  exercises  indicate  that  the  analog  approach 
to  crack-growth  prediction  possesses  some  advantages  and  some  disadvantages 
with  respect  to  digital  simulation.  Generally,  the  analog  can  accommodate 
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random-process  driving  functions  more  conveniently,  while  the  digital 
computer  is  superior  for  the  simulation  of  details. 

The  Forman  equation  was  easily  implemented  on  the  analog  for  the  special 

case  of  a growth-rate  exponent  N = 3 (representative  of  aluminum  alloys). 

F 

Formulation  in  terms  of  the  inverse  damage  parameter  **h"  makes  the  simulation 

stable  and  permits  studies  of  cracks  growing  between  0.01  and  approximately 

4 inches.  Random-process  drivers  can  be  simulated  directly  by  band-limited 

white  noise,  with  or  without  shaping  filters  to  change  the  probability 

density  function.  Simulation  of  other  materials  with  fractional  rate 

exponents  N is  possible  in  principle  by  using  the  log-antilog  components 
F 

available  on  the  analog  patchboard.  However,  these  components  are  essentially 
specialized  function-generators,  and  the  effects  of  their  possible  inaccuracies 
have  not  been  assessed  quantitatively  for  the  crack-growth  problem. 

With  regard  to  retardation  effects,  it  is  apparent  that  a digital  simula- 
tion will  be  superior.  Although  the  Wheeler  retardation  model  can  be  simulated 
on  the  analog  (again,  at  the  price  of  using  log-antilog  components) , any  attempt 
to  treat  the  more  sophisticated  load- interaction  models  forces  the  simulation 
toward  long  computing  times  because  of  the  need  for  cycle-by-cycle  detail  in 
the  stress-parameter  signals. 

Local  geometry  effects  for  two-dimensional  cracks  appear  to  be  amencible 
to  analog  simulation  by  using  a piece-wise  linear  function-generator  with  one 
independent  variable  (a  standard  patchboard  component) . A reasonable  simulation 
of  one  such  factor  was  achieved  for  a very  severe  case  involving  unbounded 
increase  of  the  stress  intensity  factor  as  a crack  grows  across  the  finite 
width  of  a plate.  Simulations  of  stable  geometry  effects,  e.g.  for  a crack 
growing  away  frcmi  a fastener  hole,  will  be  much  easier.  Again,  the  digital 
computer  is  superior  in  comparison  to  existing  Air  Force  analog  hardware  for 
the  simulation  of  geometry  factors  for  three-dimensional  cracks,  which  require 
two  input  varicibles.  However,  commercial  analog  hardware  is  available  for  this 
task  (see  Subsection  3.3). 

Analog  simulation  of  the  effects  of  GAG  cycles  has  been  shown  to  be  feasible. 
Separate  simulations  were  implemented  for  the  isolated  large  stress  range  associated 
with  each  takeoff,  and  for  the  effect  of  the  change  in  mean  stress  from  taxi  to 
flight  loads. 
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Investigation  of  the  use  of  shaping  filters  indicated  that  both  Rayleigh 
cmd  Weibull  processes  could  be  produced  from  the  band- limited  Gaussian  noise 
sources  available  on  the  analog.  However,  the  results  of  the  tests  also 
indicated  some  difficulty  in  accurately  reproducing  the  desired  shape  parameters 
for  these  distributions.  Further  tests  are  required  to  assess  the  crack-size 
inaccuracies  which  might  result  from  inaccuracy  of  the  process  shape  parameter. 

The  threshold  stress  intensity  effect  was  simulated  successfully  with  no 
difficulty.  However,  analog  simulation  of  the  threshold  effect  cannot  be 
extended  to  include  the  "shutoff’  phenomenon,  which  may  occur  when  load- 
interaction  effects  are  properly  accounted  for. 

Finally,  it  has  been  shown  that  coupled  crack  dynamics  can  be  simulated 
on  the  analog.  However,  there  is  no  reason  why  an  equivalent  digital  simulation 
could  not  be  implemented. 

In  summary,  the  results  of  the  analog  simulations  indicate  that  the  analog 
computer  provides  about  the  same  level  of  capability  as  the  digital  computer, 
but  with  some  important  differences  in  detail.  Hence,  there  appears  to  be  no 
overwhelming  reason  to  convert  from  the  digital  to  the  analog  approach  to 
crack-growth  prediction,  especially  since  many  digital  programs  have  already 
been  developed,  and  in  view  of  the  fact  that  many  structural  engineers  are 
unfcimiliar  with  analog  programming  techniques. 

3.3  Equipment  Availability 

In  recent  years  the  usage  of  existing  analog  and  hybrid  computers  and  the 
demand  for  new  machines  has  declined  signif iccuitly . Complete  listings  of  present 
day  manufacturers  of  hybrid  computers  are  difficult,  if  not  impossible,  to 
obtain.  A request  for  information  sent  to  the  six  manufacturers  listed  in  the 
1972  edition  of  the  Computer  Yearbook  yielded  only  two  responses  (see  Appendix  C) . 
Some  or  all  of  the  companies  that  did  not  respond  may  no  longer  be  in  existence. 

Fortunately  for  the  analog/hybrid  user,  the  responding  manufacturers  have 
continued  to  improve  and  refine  their  products.  In  particular  Electronics 
Associates,  Inc.  has  added  multivariable  function-generators  to  their  line  of 
products  and  accessories.  Although  digital  multivariable  function  generators 
have  been  in  use  for  some  time,  the  EAI  function-generator  employs  analog 
computation  between  function  break  points.  The  multivariable  function  generator 
acts  independently  of  the  digital  central  processor  and  thus  achieves  a combina- 
tion of  speed  and  accuracy. 
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Appendix  D compares  the  features  of  available  new  machines  with  the 
existing  CI-5000/5  computer  at  Wright-Patterson  Air  Force  Base.  Although 
some  of  the  conveniences  and  more  advanced  design  philosophies  of  today's 
new  machines  are  not  reflected  in  this  table,  the  basic  features  required 
for  solving  the  problems  discussed  in  this  report  are  listed.  The  existing 
Wright-Patterson  hybrid  computer  compares  favorably  on  all  requirements,  with 
the  exception  of  multivariable  function -generators . 
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Section  IV 


CONCLUSIONS  AND  RECOMMENDATIONS 


4,1  Summary 

The  results  of  a program  to  develop  an  efficient  approach  to  the  risk 
analysis  of  random- load  crack  growth  have  been  presented.  The  two  major 
objectives  of  the  program  were  to  implement  certain  specific  simulations  of 
crack  growth  on  hybrid  analog/digital  hardware,  and  to  develop  an  improved 
approach  to  the  modeling  of  random  loads. 

Under  the  first  objective,  eight  of  nine  specific  simulation  tasks  were 
implemented  and  verified  on  analog/digital  hardware  available  at  the  USAF 
Aeronautical  Systems  Division  computing  facility  (ASD/ADSD) . These  simula- 
tions were  achievable  well  within  the  existing  hardware  capability  (i.e., 
numbers  of  multiplier,  summation,  integrator  circuits,  etc.)  with  one  exception. 
The  exception  is  the  geometry  factor  required  for  stress-intensity  calculations 
associated  with  three-dimensional  CMCks.  Calculation  of  this  factor  requires 
a function-generator  with  two  independent  variables;  commercial  hardware  possess- 
ing this  capability  has  been  identified. 

The  ninth  task  involved  load- interaction  (retardation)  effects  and  could 
not  be  fully  implemented  as  an  analog  simulation.  The  Wheeler  retardation  model 
could  be  simulated  by  using  logarithm  and  anti— log  circuits,  but  the  accuracy 
of  these  circuits  is  questionable.  Attempts  to  formulate  the  more  refined 
cycle-by-cycle  retardation  models  were  completely  unsuccessful  because  of  the 
difficulty  in  extracting  the  discrete  information  required  for  calculating 
retardation  from  the  continuous  simulation  of  stress  time-history. 

Under  the  second  objective,  a general  procedure  has  been  proposed  for 
characterizing  random  loads  from  surveys  of  recorded  flight  data,  in  terms  of 

those  load  statistics  which  influence  crack  growth  most  significantly.  In  the 

course  of  this  development,  a general  study  of  how  to  formulate  the  problem  of 

predicting  crack  growth  due  to  random  loads  led  to  two  improvements  in  the 

approach.  First,  the  crack  growth- rate  equations  were  reformulated  in  terms 
of  a damage  parameter  inversely  related  to  crack  size.  Second,  the  "estimation 
theory" developed  by  modern  theorists  in  guidance  and  control  has  been  considered 
for  crack  growth  prediction  to  obtain  a mere  efficient  numerical  approach  to  risk 
analysis,  as  described  in  the  following  subsection. 
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4.2  Estimation  Theory  for  the  Prediction  of  Crack-Growth  Statistics 


In  the  course  of  investigating  modeling  methods  for  load  statistics,  the 
usefulness  of  estimation  theory  for  application  to  crack-growth  risk  analysis 
was  identified.  Estimation  theory  [21,22]  is  concerned  with  the  analysis  of 
the  states  of  dynamic  systems,  where  such  states  are  not  accurately  known  for 
three  reasons.  First,  the  initial  conditions  of  the  system  may  not  be  accu- 
rately known.  Second,  the  dynamic  system  is  affected  by  input  noise.  Third,  » 

measurements  of  the  state  are  also  affected  by  noise,  so  direct  observation  of 
the  state  is  not  possible.  Because  the  states  are  not  accurately  known,  they 
can  only  be  described  statistically.  A diagram  of  this  situation  is  presented 
in  the  upper  portion  of  Fig.  23.  It  should  be  noted  that  the  noises,  states, 
and  measurements  shown  in  the  figure  are  vector  quantities,  not  merely  scalars. 

It  is  convenient  for  the  purposes  of  the  estimation  theory  to  represent 
the  dynamic  system  in  "state-space"  formulation.  The  differential  equations 
of  the  dynamic  system,  of  total  order  n,  can  always  be  rewritten  as  n first- 
order  differential  equations  in  terms  of  n variables.  These  n variables  are 
collected  in  the  "state  vector"  {x}.  The  system  equation  can  then  be  written 
as 


itw 

'state 

vector" 


= {f({x})}  + [G({x})]  {w} 


"white 

noise" 


(76) 


For  the  case  of  a linear  system,  {f({x})}  is  simply  linear  in  {x},  and 
[G({x})]  is  constant,  so  that  the  system  equation  can  be  written  as 

^ {x}  = [F]{x}  + [G]{w} 


It  must  be  noted  that  the  block  diagram  of  the  system  shown  at  the  top  of 
Fig.  23  has  Gaussian  noise  input,  while  the  equations  shown  above  allow  only 
for  "white  noise"  (Gaussian  noise  with  flat  power  spectral  density) . This 
discrepancy  is  easily  handled  because  any  Gaussian  noise  can  be  produced 
by  passing  white  noise  through  an  appropriate  linear  filter.  As  shown  in 
the  lower  diagram,  an  "augmented  system"  can  be  developed,  with  an  "augmented 
state  vector"  including  the  states  necessary  to  represent  the  dynamics  of  the 
"input  noise  filter".  The  augmented  problem  is  then  appropriate  for  treatment 
using  the  equations  shown  above. 
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The  development  of  estimation  theory  began  with  consideration  of  steady- 
state  (stationary)  behavior  of  linear  systems,  with  all  variables  having 
Gaussian  statistics.  This  restricts  the  problem  represented  by  Fig.  23  in 
three  ways.  First,  the  system  equation  can  be  expressed  as  in  Eq.  77,  which 
restricts  the  system  dynamics  to  be  linear  and  the  noise  input  to  be  Gaussian. 
Second,  the  matrices  [F]  and  [G]  are  not  time-varying,  which  restricts  the 
system  and  the  noise  input  to  be  both  stationary.  Third,  the  system  has 
reached  a steady  state,  or  in  other  words  enough  time  has  passed  for  the 
initial  conditions  to  have  long  ceased  to  have  an  effect  on  the  problem.  This 
restrictive  form  of  estimation  theory  is  not  applicable  to  the  crack-growth 
problem,  which  is  an  entirely  transient  problem  without  a steady  state. 

Estimation  theory  was  then  extended  to  consider  non-steady-state  (non- 
stationary) situations,  but  the  models  of  the  system  and  noise  remained  linear 
and  Gaussian.  The  initial  conditions  of  the  system  became  important.  This  led 
to  the  development  of  the  Kalman  filter  for  estimation  of  the  state  of  a system. 
Referring  again  to  Fig.  23,  the  Kalman  filter  (measurement  filter)  keeps  track 
of  the  statistics  of  the  state  variables,  which  are  entirely  represented  by 
their  vector  of  means  and  matrix  of  covariances,  since  all  statistics  are 
assumed  Gaussian.  The  Kalman  filter  updates  these  statistics  based  on  a number 
of  effects.  First,  the  system  dynamics  cause  the  expected  value  of  the  state 
to  change,  with  corresponding  changes  in  the  covariances.  Second,  the  noise 
input  to  the  system  increases  the  uncertainty  of  the  states,  which  is  realized 
as  an  increase  in  the  covariances.  Third,  the  measurements  provide  information 
which  is  useful  in  better  estimating  the  expected  value  of  the  states.  Based 
on  the  relative  accuracy  of  the  measurements  (as  described  by  the  measurement 
covariances)  and  the  accuracy  of  the  current  estimate  of  the  state  (as  described 
by  the  state  covariances),  an  adjustment  of  the  expected  values  of  the  states 
is  made  towards  the  measured  values,  and  the  increased  accuracy  of  the  estimate 
is  realized  as  a decrease  in  the  state  covariances. 

The  Kalman  filter  has  been  extremely  valuable  in  situations  where  the  best 
estimate  of  a state  variable  was  very  important,  but  accurate  measurements  of 
that  state  variable  were  not  available.  For  instance,  navigation  of  a vehicle 
is  often  dependent  on  noisy  (inaccurate)  position  measurements,  especially  when 
operating  away  from  regular  routes.  In  such  a situation,  the  Kalman  filter  can 
greatly  improve  the  estimate  of  current  position  and  heading.  In  the  utilization 
of  the  Kalman  filter,  the  greatest  difficulty  occurs  in  accounting  for  the  measure- 
ments. The  propagation  of  the  statistics  of  the  state,  considering  the  system 
dynamics  and  the  Gaussian  noise  input,  is  straightforward. 
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Many  attempts  have  been  made  to  extend  estimation  theory  to  cases  of 
nonlinear  systems  or  linear  systems  with  non-Gaussian  statistics,  but  these 
were  not  entirely  successful.  Jazwinski  [21]  says  that 

"Using  linearization  of  one  sort  or  another,  Kalman-like 
filtering  algorithms  were  developed  and  applied  to  non- 
linear problems.  This  statistical  work  received  its 
impetus  from  the  aerospace  dollar.  Work  was  duplicated, 
triplicated;  everyone  derived  his  own  Kalman  filter, 
perhaps  partly  because  of  a lack  of  understanding  of 
Kalman's  original  work." 

It  was  soon  learned  that  development  of  estimation  procedures  for  non- 
Gaussian  statistics  and  nonlinear  systems  required  a much  deeper  understanding 
of  the  fundamental  probabilistic  structure  of  the  estimation  problem,  than  did 
the  development  of  the  Kalman  filter.  This  fundamental  probabilistic  structure 
is  now  more  completely  understood,  but  estimation  of  non-Gaussian  statistics 
with  nonlinear  dynamic  systems  is  still  not  commonly  performed  because  of  the 
immense  computational  difficulties  involved.  First,  the  means  and  covariances 
of  the  state  variables  no  longer  completely  define  their  statistics.  An  infinite 
number  of  statistical  parameters  would  theoretically  be  necessary  to  completely 
define  the  statistics.  However,  as  long  as  the  random  variables  are  defined 
so  their  statistics  are  not  far  from  normal,  a small  number  of  these  parameters 
(such  as  the  first  few  moments  of  the  probability  density  functions)  can  serve 
as  a useful  approximation  to  the  complete  statistics.  The  number  should  be 
kept  as  small  as  possible  since  each  additional  parameter  disproportionately 
increases  the  amount  of  information  that  must  be  carried  throughout  the 
computations . 

The  system  dynamics,  which  are  now  considered  to  be  nonlinear,  also 
theoretically  require  an  infinite  number  of  parameters  to  completely  specify 
their  behavior.  However,  the  forms  of  {f({x})}  and  [G({x})]  are  usually 
provided  with  sufficient  accuracy  by  closed-form  expressions.  In  the  develop- 
ment of  the  propagation  of  the  statistics  for  the  nonlinear  non-Gaussian  system, 
the  dependence  of  these  system  dynamic  descriptors  will  be  approximated  by  taking 
the  first  few  terms  of  their  Taylor  series  expansion  in  the  states,  about  the 
means  of  the  states. 
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For  the  linear  Gaussian  case,  it  was  pointed  out  that  consideration  of 
measurements  greatly  increased  the  complexity  of  the  filter.  For  the  nonlinear 
non-Gaussian  case,  this  increase  in  complexity  is  magnified  many  times  over. 
Fortunately,  treatment  of  measurements  will  not  be  part  of  the  crack-growth 
problem  being  considered. 

The  propagation  of  the  first  few  moments  of  the  probability  density  of  the 
state  for  the  nonlinear  non-Gaussian  case  is  substantially  more  complicated 
than  for  the  linear  Gaussian  case.  However,  it  is  a straightforward  calculation 
which  can  be  implemented  on  a digital  computer.  There  are  a number  of  effects, 
not  present  in  the  linear  Gaussian  case,  which  may  be  important.  First,  the 
effect  of  the  nonlinearity  of  the  system  cannot  be  approximated  by  merely 
linearizing  about  the  mean  value  of  the  states,  because  this  would  ignore  the 
very  important  variations  in  how  the  system  operates  on  other- than-average 
values  of  the  states.  Second,  the  analysis  must  consider  more  than  the  first 
two  moments  of  the  probability  density  (mean  and  covariance)  even  if  they  are 
the  only  results  desired,  because  in  the  presence  of  the  system  nonlinearities 
higher  density  moments  will  affect  the  propagation  of  the  lower  moments.  Third, 
the  non-Gaussian  statistics  of  the  initial  conditions  must  be  considered.  Fourth, 
the  non-Gaussian  statistics  of  the  input  noise  should  be  considered.  This  fourth 
effect  is  probably  the  least  important  because  the  input  noise  effects  are 
averaged  over  the  time  period  of  concern,  but  even  this  effect  should  not  be 
ignored  unless  examination  of  the  particular  problem  shows  it  to  be  unimportant. 

The  problem  of  predicting  crack-growth  statistics  can  now  be  put  in  the 
context  of  estimation  theory.  It  will  be  necessary  to  model  the  behavior  of 
crack  growth  by  a dynamic  system,  and  to  model  the  important  effects  of  the 
loading  of  the  structure  by  a random  noise  input.  It  will  also  be  necessary 
to  consider  the  dynamic  system  as  nonlinear,  and  the  statistics  of  the  noise 
and  state  variables  as  non-Gaussian. 

While  Gaussian  statistics  can  be  completely  represented  by  the  first  two 
moments  of  the  probability  density  (mean  and  covariance) , non-Gaussian  statistics 
generally  require  an  infinite  number  of  statistical  parameters  to  completely 
represent  them.  The  non-Gaussian  statistics  are  therefore  approximated  by  taking 
only  the  first  few  density  moments  (perhaps  the  first  four) . This  approximation 
will  be  reasonable  only  if  the  higher-order  moments  do  not  have  significant 
effect,  which  is  equivalent  to  saying  that  the  statistics  are  not  far  from 
Gaussian.  This  will  usually  be  the  case  if  and  only  if  the  system  dynamics 
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are  not  far  from  linear.  Unfortunately,  the  dynamics  of  crack  growth  are 
highly  nonlinear  if  they  are  formulated  with  crack  size  itself  as  a state 
variable,  as  explained  in  Subsection  2.1.  The  solution  as  presented  there 
is  to  reformulate  the  crack-growth  dynamics  in  terms  of  a damage  parameter 
"h" . As  was  shown  in  Figs.  1 and  2,  this  results  in  a linear  dynamic  system 
for  the  simple  case  of  the  Paris  crack-growth  equation,  and  a near-linear 
system  for  more  realistic  models  such  as  the  Forman  equation.  Without  the 
use  of  this  damage  parameter,  the  employment  of  the  first  few  probability 
density  moments  to  approximate  the  nonlinear  statistics  would  be  inadequate, 
and  the  use  of  a sufficient  number  of  moments  would  be  unrealistically  compli- 
cated. 

Even  with  the  necessity  to  consider  the  nonlinear  and  non-Gaussian 
situation,  the  major  difficulty  of  the  associated  estimation  theory  is  avoided 
by  the  crack-propagation  problem  being  considered  here.  This  difficulty  is 
the  treatment  of  measurements,  which  do  not  exist  in  the  prediction  problem. 
Without  measurements,  the  nonlinear  non-Gaussian  situation,  which  has  been 
either  avoided  or  implemented  only  at  extreme  cost  by  investigators  considering 
other  problems,  becomes  feasible.  An  additional  factor  which  makes  it  feasible 
to  handle  this  situation  is  that  the  computations  will  be  performed  on  large- 
scale  ground-based  computers,  while  many  applications  of  estimation  theory 
have  required  that  computations  be  performed  in  real  time  on  much  more  limited 
on-board  computers. 

The  system  dynamics  representing  the  behavior  of  crack  growth  are  relatively 
easy  to  develop.  Available  crack-growth  models  such  as  the  Forman  equation  are 
utilized,  with  the  necessary  transformation  to  model  the  behavior  in  terms  of 
a damage  parameter  "h". 

The  input  noise  to  the  dynamic  system  models  the  essential  characteristics 

of  the  loading  experienced  by  the  structure,  such  as  stress  range  As  and  peak 

stress  s . The  treatment  of  complex  effects,  such  as  retardation,  may  be 
max 

included  in  the  development  of  the  input  noise  model,  as  was  described  in 
Subsection  2.2. 

The  statistics  of  the  initial  conditions  of  the  system  (corresponding  to 
the  initial-crack  sizes)  must  also  be  modeled.  The  initial-crack  sizes  are 
usually  quite  small,  with  some  occurrence  of  larger  initial  defects.  The 
probability  density  function  of  initial-crack  size  might  very  well  appear  as 
shown  in  Fig.  24,  with  two  peaks  (a  bimodal  distribution).  This  type  of 
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density  is  not  well  represented  by  its  first  few  moments.  However,  bimodal 
distributions  can  be  handled  by  considering  the  density  to  be  the  sum  of  two 
separate  densities,  as  shown  in  the  figure.  The  two  densities  are  utilized 
as  initial  conditions  for  two  separate  analyses,  and  the  results  of  the 
separate  analyses  are  then  combined.  This  illustrates  the  capcibility  of 
handling  a multi-modal  density  function  accurately,  while  still  approximating 
the  non-Gauss ian  statistics  by  only  the  first  few  moments  of  the  density 
function. 

4.3  Discussion  and  Conclusions 

Although  eight  of  nine  analog  simulation  tasks  were  successfully  imple- 
mented, evaluation  of  the  computing  performance  indicated  that  the  analog 
approach  to  cycle-by-cycle  or  flight-by-flight  crack  growth  prediction  does 
not  achieve  any  real  improvement  in  computing  efficiency,  in  comparison  with 
equivalent  digital  computer  codes.  The  inability  of  analog  equipment  to 
process  high-frequency  signals  implies  long  computation  time  for  one  airplane 
life.  Also,  analog-circuit  errors  are  significant  enough  to  require  repetition 
of  the  analysis  several  times,  simply  to  assess  the  error  in  crack  size  for  a 
deterministic  load  history.  If  Monte  Carlo  simulation  to  obtain  the  crack- 
size  distribution  due  to  random  loading  is  attempted,  very  uneconomical 
computation  times  result,  and  the  circuit  errors  cannot  be  separated  from  the 
intended  random  effects  of  the  applied  loading.  Finally,  load-interaction 
effects  cannot  be  simulated  conveniently  with  analog  equipment.  Therefore, 
analog  Monte  Carlo  simulation  is  considered  to  be  a less  useful  approach  to 
risk  analysis  than  the  estimation- theory  formulation. 

Regardless  of  the  approach  used  for  risk  analysis,  it  has  been  shown 
that  statistical  treatment  in  terms  of  a damage  parameter  as  described  in 
Subsection  2.1,  rather  than  directly  in  terms  of  crack  size,  can  significantly 
improve  the  accuracy  of  and  reduce  the  complexity  of  the  analysis. 

The  apparent  computational  efficiency  of  estimation  theory  results  from 
its  avoidance  of  direct  numerical  summation  of  crack  growth- rate  equations. 
However,  estimation  theory  must  be  formulated  carefully  for  application  to 
the  crack-growth  prediction  problem,  for  which  the  dynamics  are  nonlinear 
and  the  statistics  are  non-Gaussian.  Computational  efficiency  will  depend  on 
a formulation  in  terms  of  near-linear  dynamics  and  near-Gaussian  statistics. 
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which  can  be  achieved  by  reformulation  of  the  growth-rate  equation  in  terms 
of  a damage  parameter.  The  growth- rate  dynamics  are  then  expanded  in  a 
Taylor  series  about  the  mean  state,  and  the  statistics  of  the  state  are 
approximated  by  the  first  few  moments  of  its  probability  density  function. 

The  numerical  procedures  for  the  application  of  estimation  theory  to 
prediction  of  crack-size  statistics  dictate  that  digital  rather  than  analog 
computation  be  adopted.  The  resulting  scheme  will  be  able  to  provide  the 
first  several  moments  of  the  probability  density  of  the  damage  parameter, 
not  merely  the  first  moment  (expected  value) . This  information  can  be  used 
to  estimate  exceedance  probabilities  of  critical  values  of  the  damage  param- 
eter (corresponding  to  critical-crack  sizes) , allowing  much  more  accurate 
risk  analysis  than  is  possible  with  present  schemes.  The  computational  costs 
are  expected  to  be  comparable  with  present  schemes  which  provide  only  the 
expected  value  of  final  crack  size. 

Some  practical  questions  about  load  models  still  remain  to  be  answered, 
i.e.  how  much  flight  data  must  be  analyzed  to  provide  a good  model  and  what 
counting  methods  should  be  used?  However,  such  questions  must  be  answered 
independent  of  the  general  approach  taken  to  the  problem  of  predicting  crack 
growth.  One  question  which  does  affect  the  utility  of  any  approach  is  how 
to  treat  load- interaction  models.  With  regard  to  the  estimation- theory 
approach,  it  has  been  shown  that  retardation  can  be  treated  efficiently,  on 
a cycle-by-cycle  basis,  by  incorporating  existing  retardation  models  in  the 
f light-data-reduction  scheme.  Appropriately  adjusted  random- load  statistics 
can  then  be  used  in  an  estimation- theory  formulation  of  the  prediction  prob- 
lem for  an  airplane  life.  Therefore,  the  estimation- theory  approach  for  risk 
analysis  of  cracks  propagating  under  random  loads  is  promising  and  warrants 
continued  development. 

4.4  Recommendations 

The  use  of  a suitably  defined  "damage  par^eter",  rather  than  crack  size 
itself,  is  highly  recommended  for  any  statistical  analysis  of  fatigue  damage, 
regardless  of  the  approach  taken.  The  slight  inconvenience  of  working  with  a 
less-than-obvious  random  variable  is  more  than  offset  by  the  improved  accuracy 
and  efficiency  obtained  with  the  damage  parameter.  Final  results  may  always 
be  transformed  for  presentation  in  terms  of  crack  size  itself,  if  desired. 
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Direct  inspection  of  large  quantities  of  actual  flight  data  is  recommended 
to  obtain  the  important  statistics  of  loading  sequences  for  utilization  in  risk 
analyses.  Although  this  entails  a substantial  effort,  the  effort  required  to 
obtain  such  statistics  when  given  an  analytically  describable  approximate  load 
model  may  not  be  significantly  less.  In  addition,  direct  inspection  of  the 
flight  data  avoids  the  errors  associated  with  analytically  describable 
approximate  models. 

Further  development  of  the  application  of  estimation  theory  to  aircraft 
fatigue  damage  risk  analysis  is  recommended.  The  techniques  of  estima- 
tion theory,  when  applied  to  analysis  of  damage  in  terms  of  an  appropriately 
chosen  damage  parameter,  promise  to  provide  efficiency  and  accuracy  in  risk 
analysis  greater  than  can  be  obtained  by  present  methods. 
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SATURATION 


FIG.  1 COMPARISON  OF  TIME-HISTORIES  OF  CRACK  SIZE  AND  DAMAGE  PARAMETER 
FOR  PARIS  EQUATION 
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SATURATION 


FIG.  2 COMPARISON  OF  TIME-HISTORIES  OF  CRACK  SIZE  AND  DAMAGE  PARAMETER 
FOR  FORMAN  EQUATION 
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FIG.  3 COMPARISON  OF  TIME-HISTORIES  OF  DAMAGE  PARAMETER  FOR  PARIS  EQUATION,  FOR 

N > 2 AND  N < 2 
P P 
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FIG.  4 


COMPARISON  OF  TIME  HISTORIES  OF  DAMAGE  PARAMETER  h = In (a)  FOR  PARIS 
EQUATION  WITH  N^  NEAR  TWO 
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FIG.  5 PROBABILITY  DENSITY  FUNCTIONS  FOR  s AND  As,  FOR  BAND-LIMITED 
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FIG.  6 REDUCTION  OF  SAMPLED  FLIGHT  DATA  TO  LOAD  STATISTICS 
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EXAMINATION  OF  SAMPLE-AND-HOLD  CIRCUITRY 
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dt 
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FIG.  8 FORMAN  EQUATION  FLOW  CHART,  USING  DAMAGE  PARAMETER  "h" 
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s and  s . 1 rms  white  noise,  s = 12.08  ksi,  s . - 0.713  ksi,  h - 10 
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FIG.  9 ANALOG  SIMULATION  OF  FORMAN  EQUATION  (REAL  TIME) 
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FIG.  10  ANALOG  SIMULATION  OF  FORMAN  EQUATION  (REPETITIVE  OPERATION) 
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FIG.  11  FINITE  PLATE  CORRELATION  FACTOR  FOR  PLATE  WIDTH  b=10  INCHES 
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FIG.  12  ANALOG  DIAGRAM  FOR  RANDOM  APPEARANCE  OF  A DETERMINISTIC 
GAG  CYCLE 
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FIG.  13  EFFECT  OF  SWITCHING  VOLTAGE  BIAS  ON  AVERAGE  TIME 
BETWEEN  GAG  CYCLES 
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FIG.  14  ANALOG  DIAGRAM  FOR  GAG/FLIGHT  LOAD  SUPERPOSITION 
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FIG.  15  SIMULATION  OF  HYPOTHETICAL  FLIGHT  WITH  THREE  SEGMENTS 
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FIG.  16  GAIN  FUNCTION  FOR  WEIBULL  DISTRIBUTION  SHAPING 
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FIG.  17  SHAPING  FILTER  FOR  RAYLEIGH  NOISE 
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FIG.  18  STRESS  INTENSITY  THRESHOLD  FLOW  CHART 
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FIG.  19  ANALOG  SIMULATION  OF  FORMAN  EQUATION  WITH  STRESS  INTENSITY  THRESHOLD 
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PIG.  20  FLOW  CHART  FOR  SIMULATION  OP  COUPLED  CRACK  GI 


FIG.  21  ANALOG  SIMULATION  OF  THE  GROWTH  OF  TWO  COUPLED  CRACKS 
WITH  DIFFERENT  INITIAL  SIZES 
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FIG.  22  ANALOG  SIMULATION  OF  THE  GROWTH  OF  TWO  INITIALLY 
EQUAL  COUPLED  CRACKS  AT  DIFFERENT  RATES 
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FIG.  23  SCHEMATIC  DIAGRAM  OF  ESTIMATION  THEORY 
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FIG.  24  TREATMENT  OF  BIMODAL  CRACK  POPULATION 
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APPENDIX  A 


FORTRAN  PROGRAM  FOR  CALCULATION  OF  CRACK-GROWTH  STATISTICS 

DTHENSIOK  DSPU(IOI)  ,DSMC  (101)  ,054(10(26)  ,GM0(26)  ,GM0E(26) 

1 READ  (5,901)  NN,nM,HAXP 

901  PCRHAT  («I10) 

IF  (NN.EO.O)  STOP 

2 READ  (5,902)  DSME AN , DSSIG1 , AZ FRO , CP 

902  FORMAT  (RFIO.O) 

WRITE  (6,911)  NN,MM,HAXP,DSMEAN,DSSIGM, AZERO,CP 
911  FORMAT  (•0NN=',I5,5X,*MM  = ' ,I5,5X,*MAXP-=' ,I5,5X,/, 

1 ' DSMEAN=' ,F15.  5,5X, 'DSSIGM=' ,F15. 5,5X,' AZERC  = ' , 

2 F15.5,5X, '0?=' ,E16.8) 

FIND  DSMO,  WHERE  DSMO  (I)  = (I  + DTH  CENTRAL  MOMENT  CF  DS 

MAXK  = FAXF  ♦ 4 
MAXKP  = MAXK  + 1 

MAXPF  = MAXP  ♦ 1 

DSMU(I)  = 1. 

DC  100  I = 2, MAXK, 2 
100  DSHO(I)  = 0. 

DSS2  = DSSIGM**2 

DS5P  = 1. 

AMOLT  = 1- 

DC  110  I = 3, MAXKP, 2 

DSSP  = DSSP*ESS2 

AMOLT  = AMOLT  ♦ FLOAT  (1-2) 

110  DSMO(I)  = AMOLT  ♦ DSSP 

WRITE  (6,913)  (DSMO(I)  ,1=1, MAXKP) 

913  FORMAT  (•  DSMO ' , / , 6 (3 X, E 1 6. 9) ) 

C 

C FIND  DSMC  WHEPF.  DSKC(I)  = (1  + 1)  TH  MOMENT  ABOUT  ZERO  CF  DS 

C 

DSMC(1)  = 1. 

FACTK  = 1. 

DC  200  K = 1,MAXK 
FACTK  = FACTK  * FLOAT  (K) 

FACTKR  = FACTK 
IR  = 0 

DSM0(K*1)  = 0. 

FACT?  = 1. 

ESMP  = 1. 

C ADD  CONTRIPOTICN  FOR  THIS  IB 

150  DSMO(K  + ‘»)  = DSMC(K*1)  ♦ DSMR  ♦ DSMU(K-IR+1) 

1 ♦ FACTK  / (FACTKR*FACTP) 

C INCREASE  IB 

IF.  = IP  + 1 

IF  (IR.GT.K)  GC  ir  170 
C PREPARE  FOR  NEW  IP  CONTP.IBOTICN 
FACTE  = FACTF  ♦ FLOAT  (IR) 

FACTKR  = FACTKR  / FLO  AT  ( K -I R ♦ 1 ) 

DSMR  = DSMF  * DSHEAN 
GO  TO  150 
170  CCNTINOE 
200  CONTINOE 
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CS4f1EN  = DSf1C(5) 

WFITE  (6,914)  (DSMO  (I)  , 1=  1 , MAXKP) 

914  FCEMAT  ('  DS MC ' , / , 6 ( 3X , E 1 6. 8 ) ) 

FIND  DS4WU  WHERE  DS4HU(I)  = (1+1) TH  CENTRAL  MOMENT  OF  (DS**4) 

DS4MU(1)  = 1. 

FACTK  = 1, 

DO  300  K = 1.MAXP 
FACTK  = FACTK  * FLOAT  (K) 

FACTKP.  = FACTK 
IF  = 0 

DS4MU(K+1)  = 0. 

FACT?  = 1. 
ns  4 ME  = 1. 

AMIF  = 1. 

C ADD  CCNTRinUTICN  FOR  THIS  IR 

250  DS4MH(K+1)  = DS4MU(K  + 1)  + AMIR  * D.S4MR  * DSMC  (4»  (K-IR)  + 1 ) 

1 * FACTK  / (FACTR*FACTKB) 

C INCREASE  TR 

IP  = IP  + 1 

IF  (IR.GT.K)  GC  TC  270 
C PREPARE  FOR  NEK  IR  CCNTPIBUTICN 
FACT?  = FACT?  * FLOAT  (IE) 

FACTKR  = FACTKE  / FLO  A*^  (K-I R+ 1 ) 

0541?  = DS4MP  * DS4MFN 
AMIS  = -AMIR 
GO  TC  250 
270  CONTINiJE 
300  CCNTIN'JE 

WRITE  (6,415)  (DS4MH  (I)  ,1=1 ,MAXPP) 

415  FORMAT  (•  ES4MIJ' ,/,6  (3X,  Elb.e)  ) 

C ’^IND  GMU  WHERE  GMU(T)  = (I  + 1)TH  CENTRAL  MOMENT  OF  G 
C 

GMEAN  = DS4MEN 
GMHd)  = 1. 
no  3000  NP  = 1,MAXP 
C 

C GIVEN  NN  = NUMBER  OF  ITEMS  BEING  SUMMED 

C NF  = ORDER  OF  MOMENT  OF  PROBABILITY  DENSITY  OF 

C SUM  WHICH  IS  BEING  COMPUTED 

C DS4Mn  = CENTRAL  MOMENTS  ITEMS  BEING  SUMMED 

C 

C DIMENSION  NFACH  .GE.  (NF/2) 

DIMENSION  NEACH(50) 

C AM  IS  THE  MOMENT  OF  THE  SUM  WHICH  IS  PEING  COMPUTED 
AM  = 0. 

C PPEPAFE  FACTNF 
PACTNP  = 1. 

DC  1090  I = 1,NP 
1090  FAOTNP  = FACTNP  ♦ FLOAT  (I) 

C 
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C FfiEPAEE  FISST  COMBINATION 
C 

K = Nr/2 

IF  (NN,LT.K)  K = NN 
IF  (K. EC-0)  GO  TO  2000 
LEVEL  = 1 
DC  1100  I = 1,K 
1100  NEACH(I)  = 2 

NEACH  (1)  = NP  - 2*  (K-1) 

C 

C PREPARE  FACTM  AND  FACTK  FOR  NEW  K 
C 

1150  FACTN  = 1. 

DO  1160  I = 1,K 

1160  FACTN  = FACTN  ♦ FLOAT (NN+ 1-1) 

C 

C COMPHTE  CONTPIBOTION  OF  THIS  CC NFIGUR ATION 
C 

1200  FACTR  = 1. 

FACTNB  = 1. 

INB  = 1 

DO  1220  I = 1,K 
KP  = NEACH  (I) 

IP  (I.F0.1)  GC  TO  1210 
INB  = INB  ♦ 1 

IF  (NEACH  (I)  . NE. NEACH  (1-1)  ) INR=  1 
1210  FACTNB  = FACTNB  * FLOAT  (INB) 

DC  1220  J = 2, KB 
1220  PACTB  = FACTB  * FLOAT  (.1) 

AMP  = FACTN  * FACTNP/  ( FACTB  ♦ FACTNB  ) 

DO  1240  I = 1,K 

1240  AMP  = AMP  * DS4H0  (NEACH  (I) +1) 

AM  = AH  + AMP 
C 

C LOOK  FOR  NEXT  EXAMPLE 
C 

IF  (K.E0.1)  GC  TO  2000 
LEVTRY  = LEVEL  + 1 
IF  (LEVTRY. EO.K)  LEVTRY  = LEVEL 
1260  IF  ((NEACH  (LEVTRY)-1).LE,N£ACH(LEVTPY*1) ) GO  TO  1400 
LEVEL  = LEVTPY 

NEACH(LEVEI)  = NEACH  (LEVEL)  - 1 
NEACH  (LEVEL *1)  = N E ACH  (I E VEL ♦ 1 ) ♦ 1 

GO  TO  1200 
C 

C REDUCE  LEVTRY 
C 

1400  LIP  = lEVTBY  + 1 
NFXCES  = 0 
DC  1270  I = LTP,K 
NEXCFS  = NEXCES  -r  NFACH(I)  -2 
1270  NEACH(T)  = 2 

NEACH  (LEVTRY)  = N E ACH  (LE VTP,  Y)  + NEXCES 
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LEVTPY  = LEVTFY  - 1 
IF  (LEVT'’Y-GE*  1)  GC  TC  1260 

C 

C LEVTRY  REACHED  ZERO  

C WHICH  MEANS  IT  IS  TIKE  TC  REDOCE  K 

K = K 1 

TF  (K.EO.O)  GC  TO  2000 
NEACH(I)  = REACH  (1)  + 2 
lEVEL  = 1 
GC  TO  1150 


C 

2000  CONTINOE 

3000  GMU(NP*1)  -•  am  / FLO  AT  (NN)  **N  P 

WRITE  (6,916)  (GMU  (I)  ,I=1.MAXPP) 

916  FORMAT  ('  GHU * ,/, 6 (3 X , El 6. 8) ) 

C FIND  GHtlE  WHICH  IS  THE  'ECCENTRICITIES'  OF  G 
C GrlUE(I)  = (T+1)TH  ECCENTRICITY 

C 


3500 

917 


922 


GMOE(I)  = 0. 

GM0E(2)  = 0. 

GH0E(3)  = C. 

AMULT  = 1. 

RGM02  = (GMU  (3)  **-  5) 

CO  3500  K = 3,KAXP,2 

GM0E(K+1)  = GMU  (K+1)  / (rGMU2**K) 

IF  ( (K  + 1) .GT.MAXP)  GC  TC  3500 
AMULT  = AMULT  ♦ FLOAT (K) 

GMUE(K  + 2)  = GM'J(K  + 2)  / (RGMU2*»  (K+ 1)  ) 

CCNTINtiE 

WPITE  (6,917)  (GMUE(I) ,I=1.MAXPP) 
FORMAT  (•  GMUE* ,/,6 (3X,E16.8)  ) 

DHMEAN  = -CF  * MM  ♦ GMEAN 
DHMn2  = (CP*FLCAT  (MM) ) **2  * GMU  (3) 
WRITE  (6,922)  GM E AN , CHM E AN , DHMU2 
FORMAT  ('  GMEAN=' ,F15.5,5X,'DHMSAN=' , 


- AMULT 


F15.5,5X,' DHMU2= 


GC  TC  1 
END 


,F15.5) 
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APPENDIX  C 


COMMERCIAL  HYBRID  COMPUTER  MANUFACTURERS 


MANUFACTURER 

COMMENT 

ADAGE,  INC. 

1079  Commonwealth  Ave. 

Brighton,  Mass.  02159 

Mostly  Display  Equipment 

APPLIED  DYNAMICS  COMPUTER  SYSTEMS 

2275  Pratt  Road 

Ann  Arbor,  Mich.  40104 

DELTA  FOUR 

Plug  & Logic  to  IBM/ 360  & XDS/SIGMA 

Optimization  & Simulation 

COMDYNA  INC. 

N.A. 

Model  648 

Plug  & Logic  to  EAI 

Simulation  & Student  Training 

ELECTRONIC  ASSOCIATES,  INC.  (EAI)* 

185  Normouth  Parkway 

West  Long  Branch,  N.J.  07764 

590,  690,  7945 

Simulation 

GPS  INSTRUMENT  COMPANY,  INC. 

188  Needham  Street 

Newton,  Mass.  02164 

200T  System 

SYSTRON-DONNER  CONCORD  INSTRUMENT  DIVISION* 

888  Galindo  Street 

Concord,  California  94520 

Model  SOY 

Simulation  & On-Line  Control 

* 

Manufacturers  which  responded  to  survey. 
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APPENDIX  D 


HYBRID  COMPUTER  FEATURES 


Item 

EAI  Pacer 

500-700  Series 

Sy s tron-Donner 

System  80 

Wright-Patterson 
Astrodata  Hybrid 
with  CI-5000/5 
Analog 

Integration/Summation 

Units 

8-42 

10-24 

85 

Summation  Units 

8-54 

36-112^^^ 

85 

Multiplier  Units 

2-60 

4-114 

171 

Single  Variable 

Function  Generators 

0-36 

3-15 

37 

Comparators/ 

Switching  Units 

2-36 

9-24 

56 

Analog  Multi-Variable 
Function  Generators 
Available 

Yes<^> 

No 

No 

Digital  Core 

8K-32K 

? 

Maximum 

Unknown 

Comments 

Wide  variety 
of  peripherals 
and  accessories 
available 

Includes 

UNIBUS  PDP-11 

central 

processor 

(1)  Numbers  dependent  upon  specific  model  and/or  series  selected.  Range 
of  values  given. 

(2)  Number  of  operational  amplifiers  given.  May  be  used  as  integrators 
or  summers. 

(3)  Includes  inverters  and  resolvers. 

(4)  Some  computation  performed  digitally. 
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